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Preface 


Although  x-ray  air  build-up  factors  can  be  an  extremely  Important 
tool  In  simplifying  the  x-ray  air  transport  problem,  no  complete  set  of 
bulld-up  factors  exists  In  current  literature.  Many  computer  codes  also 
exist  to  handle  the  x-ray  air  transport  problem,  but  none  have  been  used 
to  calculate  a complete  set  of  bulld-up  factors. 

I have  used  the  least  expensive  code  which  retains  adequate  accuracy 
to  compile  a complete  sot  of  x-ray  air  bulld-up  factors.  I must  acknowledge 
the  heir  of  Dr.  C.  J.  Bridgman  for  his  help  In  guiding  me  In  the  correct 
direction  for  making  this  study  a success.  I also  must  acknowledge  the 
help  of  Major  Winfield  S.  Bigelow  who  was  Invaluable  In  getting  the  program 
to  run  successfully. 
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AFIT/CKE/FH/?6l>-7 


This  report  is  a compilation  of  time  integrated  x-ray  energy  build-up 


factors  from  a monoenergetic  point  source  in  infinite  homogeneous  air 


These  factors  were  computed  by  the  use  of  FH0TDI5,  a moments  method  computer 


code,  and  performed  on  a CDC  6600.  This  code  was  chosen  ufter  a literature 


search  and  a review  of  many  computer  models.  Energies  from  12  keV  to  1000 


koV  and  ranges  from  1 mean-free-path  to  15  mean- free- paths  are  considered 


All  results  ara  presented  on  semi-log  graphs  with  each  graph  containing 


one  energy.  This  program  is  estimated  to  have  an  error  of  at  most  20#. 

The  results  are  compared  to  Eonte  Carlo  and  Discrete  Ordinates  calculations 


Even  though  the  moments  calculations  do  not  completely  agree  with  any  of 
the  Eonte  Carlo  calculations,  the  moments  calculation  agree  with  the  average 


of  the  Eonte  Carlo  calculations.  A complete  derivation  of  the  moments 


method  from  the  Boltzmann  Transport  Equation  is  also  included 


X-RAY  BUIID-UP  FACTvliS 


I iEtiodu.cUfln 


Purpose 

The  purpose  of  this  report  is  to  provide  one  complete  set  of  time 
integrated  monoenergetic  x-ray  energy'  build-up  factors  in  infinite  horcoge- 
eous  air.  This  study  has  used  the  moments  method  in  the  form  of  the 
program  PHOTDIS  (Ref  l).  This  program  is  written  in  Fortran  and  was  run 
on  the  CDG  6600  of  the  Aeronautical  Systems  Division  at  Wright-Patterson 
Air  Force  Base,  Ohio.  This  study  was  undertaken  because  of  the  lack  of  a 
complete  set  of  these  build-up  factors. 

Jhsrcx  of  £y,Ur.;-us  Jfcslaxa 

The  normal  treatment  of  photon  attenuation  in  an  absorbing  medium  at 
range  r from  an  isotropic  monoenergetic  point  source  is 

- -U* 

F(r)  ■=  i55— - (1) 

4trr 

where 

F(r)  is  the  energy  fluence  at  distance  r in  Joules/ra  (or  Cal/cni  ) 

S is  the  total  energy  emitted  by  the  monoenergetlc  source  in  Joules 
(or  Calories) 

u is  the  total  macroscopic  cross  section  in  m ^(or  cm  *) 
r is  the  distance  of  interest  in  meters  (or  centimeters) 

This  equation  accurately  describes  the  fluence  if  the  x-rays  interact  only 
by  photoelectric  absorption.  So  this  equation  is  accurate  for  low  energy 
x-rays,  energy  below  about  12  keV  where  photoelectric  absorption  dominate. 


i 


% 


For  x-rajs  of  energy  above  12  koV,  Conpton  scatter  play3  In  increasingly 
important  part.  Thoreforo,  Eq  (l)  must  be  modified  to  account  for  the 
Compton  scatter.  One  possible  modification  is 

pc  "Mr 

F(r)  - ■ES2-t-  (2) 

4nr2 

where  B i3  the  energy  dependent  build-up  factor,  or  the  correction  factor 
to  account  for  Compton  scatter. 

Eq  (2)  assumes  a monoencrgetic  source  which  is  not  realistic  in 
nuclear  effects  calculations.  However  the  real  polycnorgetic  sources 
can  be  treated  by  using  Eq  (2)  in  a multienergy  group  calculation! 


(3) 


where  represents  the  fraction  of  the  continuous  source  energy  spectra 
emitted  within  tho  bounds  of  energy  group  i.  Eq  (2)  also  assumes  that 
the  x-ray  fluence  varies  only  with  one  spatial  cooixUnato  which  becomes 
increasingly  Invalid  at  larger  energies  and  higher  altitudes.  This  assump- 
tion is  made  to  prevent  the  calculational  complication  of  a two  dimensional 
geometry.  However,  two  dimensional  effects  can  be  approximated  by  rearranging 

o 

Eq  (3)  into  a 4rrr  fluence  and  employing  mass  Integral  scaling  (Ref  2)1 

%)T2F  » ^ ^ exp(-/  iijT'dr')  (4) 

These  equations  also  Ignore  the  tine  variable  of  the  x-ray  because 
most  targets  respond  to  total  x-ray  dose  not  x-ray  dose  rate.  Finally,  tho 
calculations  for  build-up  factors  are  presented  in  terns  of  mcan-free-path 
which  makes  them  applicable  to  homogeneous  air  at  any  altitude. 

■SaflKC 

Duo  to  the  limits  of  tho  convergence  of  tho  moments  method  program,  the 


x 


ranges  considered  were  from  one  to  15  mean-free-paths  of  the  source  photon 
energy*  The  build-up  factor  for  energies  below  about  12  keV  is  nearly  one, 
so  energies  below  12  keV  were  not  investigated.  This  program  does  not 
include  pair  production,  so  energies  used  in  the  calculations  are  all 
below  1000  keV.  To  facilitate  the  users'  evaluation  of  the  accuracy  of 
the  build-up  factors  presented  here,  comparison  to  other  calculations  are 
also  included  in  this  report. 

Plan  of  Development 

This  report  starts  with  an  explanation  of  computer  models  which  are 
used  in  x-ray  transport  problems.  A discussion  of  the  program  PHOTDIS 
follows.  The  graphs  of  the  build-up  factors  for  each  energy  are  the  main 
portion  of  this  report.  The  discussion  of  the  results  contain  a comparison 
to  others  who  have  made  x-ray  transport  calculations.  Finally  a detailed 
description  of  tho  derivation  of  the  moments  equation  and  the  reconstruction 
of  the  energy  fluence  can  bo  found  in  the  appendix. 
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II 


Three  numerical  models  Jiave  evolved  which  are  applicable  to  compute 
total  fluence  of  x-rays  where  energy  is  greater  than  12  keV.  The  three 
models  are  Monte  Carlo,  Discrete  Ordinates,  and  the  Moments  Method.  Each 
modol  has  advantages  and  disadvantages  which  ’..’ill  be  discussed  along  with 
a brief  description  of  each  method.  Since  the  moments  method  was  selected 
as  the  model  used  in  this  project,  it  will  be  presented  in  depth  in 
Appendix  A as  well  as  the  briof  overview  in  this  section. 


In  the  Monte  Carlo  method,  a photon  leaving  the  source  is  traced  through 
the  medium  of  interest.  As  the  photon  is  traced,  a record  is  kept  of  its 
energy,  direction  and  position.  Compton,  pair  production  and  photoelectric 
events  are  experienced  by  the  photons  in  accordance  with  the  cross  sections 
for  these  events  at  a rate  determined  by  statistical  probability.  The 
photon  is  followed  until  it  is  absorbed  or  until  it  passes  out  of  the  area 
of  interest.  Many  photons  are  traced  by  this  method.  At  a predetermined 
distance  r,  the  energy  of  all  the  photons  reaching  that  distance  is  added 

p 

together  and  divided  by  4nr  . Thie  is  the  total  energy  fluence  at  that 
distance.  To  minimize  computations  without  sacrificing  accuracy, 
statistical  Improving  achemes  are  used.  These  schemes  aro  used  to  sample 
only  "important"  photons,  which  are  determined  by  a "weight".  The  weight 
of  a photon  is  an  artifical  biasing  which  is  introduced  while  tracing  the 
photon.  Several  statistical  improving  schemes  axe  statistical  estimationsi 
exponential  transformations,  Russian  Roulette  and  non-nbsorptivo  weighting. 
(Ref  3119-20) 

Statistical  estimation  may  be  a lost-flight  estimator  which  calculates 
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tho  probability  that  a photon  which  hue  Just  undergone  a Compton  event 
will  travel  a predetermined  distance  bofore  absorption.  Another  otntlotlcal 
estimation  nchomo  la  noxt-ovout-estimator.  This  scheme  calculates  the 
probability  that  a photon  which  has  Just  scattered  will  roach  a predetermine <1 
dlatance  after  It  la  involved  in  another  oc&ttor.  (Hof  3*20) 

Exponential  transformation  Mason  tho  distance  to  tho  next  evont  to 
increase  tho  probability  of  n long  flight.  Tho  wolght  associated  with  a 
photon  la  adjusted  to  correct  for  the  bins.  (Ref  3*21 ) 

If  a photon  lun  lost  no  much  of  Its  weight  through  adjustment  for 
biasing*  Russian  Roulette  may  bo  used.  A random  number  Is  generated  and 
compared  with  a survival  probability,  which  is  also  less  than  one.  If  the 
random  number,  which  1s  also  less  than  ono,  Is  I033  than  tho  survival 
probability,  tho  weight  la  multiplied  by  tho  reciprocal  of  tho  survival 
probability  anl  the  tracking  is  continued.  If  tho  random  number  lu  greater 
than  tho  survival  probability,  tho  tracking  is  discontinued.  (Ref  3*21 ) 

In  non-nbsorptivo  weighting,  at  each  event  the  photon's  weight  la 
sciilcd  by  a probability  tluxt  the  photon  will  not  be  absorbed.  Tho  scaling 
factor  is  tho  scattering  croon  section  divided  by  the  total  cross  noction. 
When  weight,  correction  reduces  the  photon's  ability  to  contribute  to  a very 
small  amount,  Russian  Roulette  is  used.  (Ref  3*22) 

The  accuracy  of  the  answer  depends  on  many  factors.  The  statistical 
fluctuation  produced  by  tho  stochastic  txaturo  of  the  method  along  with  the 
typo  of  weighting  used  nro  two  major  factors.  Other  factors  affecting  tho 
accuracy  nro  the  cross  sections  and  tho  number  of  photon  tracked  (histories). 
Tho  number  of  histories  Is  limited  by  the  slsr.o  of  tho  computer  ami  the 
amount  of  computer  time  and  money  spent,  on  the  problem.  Rut  the  Monte  Carlo 
method  can  treat  a complex  goomotry  consisting  of  many  different  materials. 
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Discrete  Ordinates  Method 


The  Discrete  Ordinates  method  deals  with  a numerical  solution  to  the 
Boltzmann  transport  equation.  The  Boltzmann  transport  equation  (which  also 
describes  time  Independent  fluence)  Is 

ft-VF(r^.E)  ♦ pMr.O.E)  - S(r,Q,E)  ♦/Jn,F(rfn,,E')  dO*  dE*  (5) 

where 

F Is  the  fluence 

fl  is  the  direction  vector 

r Is  the  coordinates 

E is  the  energy 

S Is  the  source  function 

Is  the  total  cross  section 

Ms  is  the  scatter  cross  section  from  direction  O'  and  energy  E'  to 
direction  ft  and  energy  E 

The  most  popular  Discrete  Ordinates  programs  solve  the  Boltzmann 
transport  equation  numerically  by  a method  known  as  the  method.  The 

method  expands  the  scatter  cross  section  In  the  Integral  term  in  Legendre 
polynomials.  By  using  the  addition  theorem  for  the  Legendre  polynomial, 
that  integral  reduces  to  one  in  energy  and  one  in  angle.  Numerical  quadrature 
is  employed  to  evaluate  the  angle  Integral.  The  equation  is  separated  into 
a group  of  equations  by  replacing  tha  energy  Integral  with  a summation  of 
group  to  group  scatter  terms.  Since  the  resulting  equations  are  differential 
equations,  finite  difference  approximations  are  used  to  redace  these  equations 
to  a group  of  coupled  algebraic  equations.  The  numerical  quadrature  and 
finite  difference  approximations  Introduce  truncation  error.  On  the  other 
hand  if  more  terms  in  the  quadrature  and  a finer  mesh  are  used,  round-off 
error  Is  increased  and  the  answer  may  not  converge.  So  some  error  is 
unavoidable.  The  computation  takes  less  time  than  Monte  Carlo  computations 
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but  can  not  handle  complex  geometries. 
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Koncnta.  Method 

The  Moments  method  also  produces  a solution  to  the  Boltzmann  Transport 
Equation  by  numerical  methods  but  Involves  more  analytical  reduction  than 
does  Discrete  Crdlnates.  In  this  method,  the  scatter  cross  section  and  the 
fluence  in  all  terms  of  Eq  (5)  are  expanded  In  Legendre  polynomials.  Using 
the  addition  theorem,  the  scatter  multiple  Integral  reduces  to  a single 
Integral  In  energy.  By  multiplying  through  the  equation  by  the  appropriate 
powerset  directional  angle  factors  and  Integrating  each  term  over  direction, 
the  equation  Is  reduced  to  a diffcrentlal-integro  equation  set.  By 
multiplying  this  set  of  equations  by  the  spatial  moments  and  Integrating 
over  all  space,  the  equations  reduce  to  the  moments  equations.  The  moments 
equations  are  rocurclvo  and  because  of  the  Integration  in  energy,  they  must 
be  evaluated  by  numerical  quadrature  in  energy  space.  Onco  the  moments 
have  been  calculated,  the  flucnce  Is  reconstructed  by  means  of  summing  all 
the  moments  of  a particular  energy  after  they  have  been  multlplyed  by  a set 
of  biorthogonal  polymonials.  A complete  derivation  of  the  moments  equation 
and  the  reconstruction  of  the  fluence  Is  given  in  Appendix  A.  (Ref  liJ-14) 

Since  only  a finite  number  of  moments  can  be  used,  this  method  is  also 
limited  In  accuracy.  But  depending  on  how  many  moments  are  used,  the  error 
can  be  limited.  Thus  a predetermined  accuracy  can  be  reached.  This  method 
Is  also  limited  to  the  time  independent  system  and  to  e simple  geometry  and 
a homogeneous  infinite  medium.  The  time  used  to  evaluate  the  fluence  by  this 
method  is  usually  less  than  either  Konto  Carlo  or  Discrete  Ordinates  in  such 
simple  geometries. 


■» 
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Ill  Program  Used 
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Lltsiaturs  Search 

Monte  Carlo.  Many  Monte  Carlo  programs  were  considered  for  thl3 
compilation  of  x-ray  build-up  factors.  The  MASTER  program  file  (Ref  4)  is 
a set  of  programs  and  cross  section  libraries  used  to  solve  radiation 
transport  problems.  The  programs  in  this  file  are  FASTER,  BETA  and  TEMPER. 

All  the  programs  can  handle  complex  geometry,  nonlinear  and  time  dependent 
photon  five nee s.  But  the  generality  of  this  file  makes  it  difficult  for 
the  user  not  acquainted  with  this  file  to  provide  data  for  a problem.  The 
DART  ar.d  DART  II  codes  (Ref  5 &rd  6)  are  codes  which  handle  time  dependent 
photon  transport  in  air.  DART  is  a one  dimensional  code  which  assumes 
homogeneous  ail,  while  DART  II  assumes  nonhomogeneous  air  and  takes  into 
account  the  curvature  of  the  earth.  DART  II  is  esplcally  suited  for  high 
altitude  air  transport.  The  THISTLE  code  (Ref  7)  is  a time  dependent  code 
which  describes  x-ray  transport  in  exponential  atmosphere. 

The  MORSE  code  (Ref  8)  is  a gamma  ray  transport  code  which  use3 
multigroup  cross  sections.  This  code  can  handle  three  or  one  dimensional 
problems.  It  treats  the  atmosphere  as  a homogeneous  medium  and  the  answers 
are  in  a time  dependent  form.  The  HAM  code  (Ref  9)  Is  a modified  MORSE  code 
which  incorporates  varying  air  density  at  high  altitudes  and  to  take  into 
account  the  curvature  of  the  earth.  One  of  the  most  popular  Monte  Carlo 
codes  used  for  x-ray  radiation  transport  is  the  PHOTRAN  code  (Ref  10,11  and  12). 
This  code  handles  time  dependent  radiation  transport  and  calculates  energy 
deposition,  photon  flux,  electron  flux  or  tissue  equivalent  dose.  It  can 
consider  energies  from  zero  to  100  MeV  taking  into  account  coherent  and 
incoherent  scattering,  photoelectric  effect,  pair  production,  fluorescent 
and  annihilation  radiation.  It  can  consider  three  dimensional  geometry  in 
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a homogeneous  medium.  ANDY  is  a series  of  transport  codes  (Ref  13)  designed 
for  time  dependent  photon  transport  which  can  handle  three  dimensional 
homogeneous  atmosphere. 

All  the  codes  listed  above  were  considered  hut  rejected  because  of  the 
amount  of  time  and  money  involved  in  Implementing  and  running  these  programs. 

Discrete  Ordinates.  DTFXKAY  is  a one  dimensional  code  (Ref  14)  using 
Discrete  Ordinates  method  of  solution.  This  code  does  not  treat  time 
dependence  and  assumes  a homogeneous  atmosphere.  ANISN  is  a Discrete 
Ordinates  code  (Ref  3)  which  treats  a homogeneous  atmosphere  in  one  dimension. 
It  dees  not  handle  time  dependence. 

Moments  Method.  PHOTDIS  i3  a one  dimensional  code  (Ref  l)  for  computing 
transport  in  a homogeneous  atmosphere.  It  does  not  liandle  time  dependence 
hut  it  does  compute  the  build-up  factors.  This  program  was  selected  because 
of  the  ease  of  input  and  the  direct  output  of  build-up  factors.  The  speed 
of  this  program  was  another  factor  in  its  selection  since  tine  dependence  and 
inhomogeneous  atmosphere  are  not  considered.  Another  factor  for  choosing  a 
moments  method  program  is  the  approximation  of  the  error  through  convergence 
analysis  which  is  not  available  from  the  other  methods  of  solution. 

Description  of  .PHOTDIS 

The  program  PHOTDIS  consists  of  two  phases*  Phase  I and  Phase  II.  In 
Phase  I,  the  moments  for  the  scattered  fluence  are  calculated.  Most  of 
Phase  I is  concerned  with  calculating  the  integral  in  Eq  (28),  which  is 
found  in  Appendix  A.  Phase  II  reconstructs  the  fluence  from  the  moments 
using  Eqs  (31). (32)  and  (33)*  The  majority  of  Phase  II  is  a numerical 
calculation  of  Eq  (33) » which  is  made  difficult  since  only  values  of  W. 
are  known.  At  the  end  of  Phase  II,  the  build-up  factor  is  calculated  using 
Eqs  (3^), (35)  and  (36).  A more  complete  description  of  the  program  can  be 
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found  in  Ref  1. 


Inrut 

The  input  to  the  program  consists  of  12  different  sets  of  data-  They 

are  i 

1)  Klein-Nlshina,  Photoelectric  and  Compton  scatter  cross  section  as  a 
function  of  energy. 

2)  Parameter  determining  if  mesh  parameters  are  to  be  input. 

3)  Altitude. 

4)  Density  of  the  air  at  the  altitude  considered. 

5)  Energy  above  which  no  correction  is  to  be  made  to  the  Klein-Nlshlna 
scattering  equation. 

6)  Parameters  determining  which  phases  are  to  be  run. 

7)  Number  of  space  points,  angular  variables  and  a parameter  determining  if 
angular  fluence  is  to  be  calculated. 

8)  Number  of  moments  to  be  used. 

9)  Source  energy  and  source  strength. 

10) Nesh  parameters. 

11) Anglos  at  which  tho  fluence  is  to  be  calculated. 

12) Distance  at  which  build-up  factors  are  to  be  calculated. 

The  cross  sections  used  for  this  study  are  from  UCRL-50174  (Ref  13)  and  from 
AFWL-TR-67-tl  (Ref  14).  A listing  of  the  code  and  a sample  input  are  furnished 
in  Appendix  B. 
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IV  Results 


Choice  of  Options 

Energies.  As  stated  in  Section  I,  the  magnitude  of  the  bulld-up  factor 
Is  related  to  the  amount  of  Compton  scatter.  For  energies  where  Compton 
scatter  is  a major  contributor  to  the  total  cross  section,  the  bulld-up 
factor  Is  greater  than  one.  Therefore,  this  study  begins  at  12  keV  since 
below  that  energy,  Compton  scatter  is  negligible.  For  energies  below  12  keV, 
Eq  (l)  can  be  used  since  the  build-up  factor  is  approximately  one.  The 
Compton  scatter  cross  section  begins  to  decrease  for  energies  above  100  keV. 
This  decrease  continues  until  the  Compton  scatter  is  very  low,  but  is  the 
only  major  contributor  to  the  total  cross  section.  Pair  production,  which 
is  zero  below  1020  keV,  increases  very  rapidly  just  above  1020  keV  and 
quickly  becomes  the  major  factor  in  the  total  cross  section.  Since  the 
program  does  not  handle  pair  production,  this  study  terminates  at  1000  keV. 

Ranges.  For  each  energy,  calculations  were  made  of  the  build-up  factor 
for  a distance  from  the  source  from  one  to  15  mean-free-paths  In  increments 
of  one  mean-free-path.  The  answers  axe  not  converged  well  beyond  15  mean- 
free-paths.  The  build-up  factors  at  the  various  distances  in  mean-free-paths 
are  valid  at  any  altitude,  but  the  mean-free-path  varies  with  altitude.  So 
at  any  altitude  of  interest,  the  build-up  factor  can  be  obtained  for  any 
distance  in  meters  by  converting  the  range  in  meters  to  mean-free-paths  at 
any  altitude.  Fig  0 shows  the  variance  of  the  mean-free-path  with  altitude 
for  different  energies. 

Grachs  of  Bulld-up  Factors 

Fig  1 to  Fig  105  are  the  graphs  of  the  build-up  factors  for  various 
energies  ranging  from  12  keV  to  1000  keV.  For  energies  between  12  and  100 
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kev,  the  increment  is  2 keV.  For  energies  between  100  and  200  keV,  the 
Increment  is  5 keV.  For  energies  between  200  and  500  keV,  the  increment  is 
10  keV  and  for  energies  higher  than  500  keV,  the  increment  is  50  keV.  On 
each  graph,  the  points  shovr.  tire  the  calculated  results  with  the  cur/e  fitted 
by  spline-fitting  technique  performed  by  the  computer  and  drawn  by  a Calcomp 
plotter  Kodel  765* 

Empirical  Build-up  Factor  Equation 

Taylor  (Ref  19)  has  developed  a simple  equation  to  describe  build-up 
factors  previously  calculated.  This  equation  is 

B-a/1*-  (6) 

where 

y is  number  of  mean-free-paths  of  source  energy 


A ,c  and  c are  constants  to  be  determined  from  calculated  Build-up  factors. 

A calculation  was  performed  using  the  build-up  factors  obtained  in 
this  study  to  obtain  these 
maximum  percent  difference 
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constants.  A list  of  these  constants  and  the 
is  shown  in  Table  I.  (Text  continues  on  page  121 ) 


Table  I. 

Constants 

for  empird  cal  build-up  factor  equation 

Energy  in  keV 

A1 

A2 

C1  °2 

Maximum  % difference 
below  8 KFP 

12 

-0.227 

1.227 

-0.400  0.000 

0.35  % 

14 

-0.370 

1.370 

-0.400  0.000 

0.89  % 

16 

-O.323 

1.323 

-0.680  0.020 

1.78  % 

18 

-O.634 

1.634 

-0.460  0.020 

0.91  % 

20 

-1.072 

2.072 

-O.36O  0.020 

2.?1  % 

22 

-1.048 

2.048 

-0.480  0.040 

2.99  % 

24 

-1.740 

2.748 

-O.340  0.040 

0.?5  % 

26 

-2.673 

3-673 

-0.260  0.040 

I.85  % 

28 

-2.664 

3-664 

-0.300  0.060 

3.87  $ 

30 

-6.038 

7.038 

-0.140  C.040 

2.90  % 

32 

-8.805 

9.805 

-0.100  0.040 

1.4  5 % 

34 

-6.504 

9.504 

-0.100  0.060 

I.63  % 

36 

-75.83 

76.83 

0.000  0.020 

0.97  % 

38 

-20.03 

21.03 

-0.020  0.060 

3.77  % 

40 

-16.94 

17.94 

-0.020  0.080 

6.93  % 

45 

14.59 

-13.59 

0.120  -0.020 

13.7  % 

50 

11.31 

-10.31 

0.160  -0.040 

21.1  % 

55 

109-2 

-108.2 

0.120  0.100 

18.0  % 

60 

-II.05 

12.05 

0.000  0.200 

26.6  % 

70 

-114.1 

115.1 

0.140  0.160 

25.7  % 

80 

-113-1 

114.1 

0.160  0.180 

27.6  % 

100 

-10.93 

11.93 

0.060  0.260 

36.O  % 

120 

-8.153 

9.153 

0.020  0.280 

39.1  % 

150 

13*14 

-12.14 

0.260  0.100 

36.3  % 

200 

-88.92 

89.92 

0.180  0.200 

33-1  % 

250 

-6.3O8 

7.308 

0.000  0.260 

39-9  % 

300 

19.89 

-18.89 

0.200  0.120 

33.1  % 

350 

72.37 

-71.37 

0.160  0.140 

31.2  % 

400 

-6.063 

7- 063 

0.000  0.220 

35-1  % 

500 

16.29 

-15.29 

0.160  0.080 

28.8  % 

600 

-57-58 

57.58 

0.100  0.120 

25.7  % 

750 

-17.20 

18.20 

0.060  0.120 

22.6  % 

14 
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V Discussion 


Intrinsic  Eatliaats  at  .tbs  Accvgasi 

To  examine  the  accuracy  of  this  program > the  convergence  of  a finite 
number  of  moments  must  be  examined  first.  If  all  the  moments  could  be 
calculated  from  Eq  (28),  found  In  Appendix  A,  no  error  from  this  source 
would  be  introduced.  To  estimate  the  error  injected  by  the  truncation  of 
terms,  examination  of  the  change  of  the  values  of  the  expansion  coefficient 
F"(r  ,E)  as  more  moments  are  used  in  Eqs  (31)  and  (32)  would  be  used.  This 
analysis  was  performed  by  Bigelow  (Ref  li 116-119)  on  this  program  and  his 
findings  indicate  that  the  total  error  produced  by  the  truncation  of  all 
but  the  first  nine  terms  results  in  a 5$  error  in  the  expansion  coefficients. 

The  error  introduced  by  calculating  the  scattered  fluence  by  using  the 
polynomial  reconstruction,  Eq  (3*0,  is  about  1$.  The  numerical  error  intro- 
duced by  the  calculations  performed  in  Eq  (l),  (29),  (35)  and  (36)  are  from 
I56  at  small  distances  for  small  energies  to  15$  at  large  distances  for  large 
energies.  The  figures  used  above  were  also  obtained  by  analysis  performed 
by  Bigelow  (Ref  li 126- 127). 

Therefore,  the  total  maximum  error  for  the  build-up  factors  introduced 
from  all  sources  is  20$,  which  only  occurs  at  large  distances  for  large 
energies. 

Comparison  of  Results  to  Others 

Introduction.  Because  of  the  limited  amount  of  data  for  build-up  factors 
of  the  energy  range  considered  in  this  6tudy  for  air,  some  modification  of 
the  data  obtained  from  other  sources  had  to  be  made.  In  seme  cases,  the 
data  obtained  from  other  sources  were  in  the  form  of  energy  fluence  instead 
of  build-up  factors.  To  obtain  the  build-up  factors  in  this  case,  Eqs  (l) 


and  (36)  were  used.  In  other  cases,  the  results  were  presented  as  a fraction 
of  vacuum  fluence  that  reached  a certain  distance.  Then,  the  following 
calculations  were  performed.  The  vacuum  fluence  can  be  defined  as 

Fv  - s/4rrr2  (7; 

By  dividing  Eq  (2)  by  Eq  (7)  the  fraction  of  the  vacuum  fluence  V is 
obtained  1 

V - Be*111  (8) 

is  equal  to  the  number  of  mean-free-paths  y.  Solving  Eq  (8)  for  B 

B - V/e~y  (9) 

So  for  any  mean-free-path  of  interest,  the  build-up  factor  can  bo  extracted 
from  data  giving  only  fraction  of  the  vacuum  fluence  to  that  distance. 

All  of  the  comparisons  made  below  are  in  terms  of  percent  difference 
which  is  defined  by  the  following  equation 1 

to  £8bbm:Uqd1  . 100J6  _ (10\ 

(Foments  method  result)  ^ u ' 

where  PD  is  the  percent  difference. 

Monte  Carlo.  Banks  (Ref  15),  using  a modified  CGRE  Monte  Carlo  code, 
performed  photon  transport  calculations  for  energies  from  20  to  300  keV. 

He  presented  some  of  the  results  in  the  form  of  a graph  of  build-up  factors 
at  10  mean-free-paths.  This  graph  as  well  as  a similar  graph  of  the  build- 
up factors  calculated  by  the  moments  method  are  superimposed  and  presented 
in  Fig  106.  The  percent  difference  for  the  energies  below  50  keV  are  as 
low  as  -47?S.  But  above  50  keV,  the  percent  differences  range  up  to  26$ 
with  the  average  being  23$.  Taking  into  account  the  error  of  the  build-up 
factor  calculation  associated  with  the  moments  method,  the  rest  of  the  error 
can  be  attributed  to  the  use  of  different  cross  sections. 

Coleman  did  a study  computing  x-ray  transport  and  the  results  are 
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presented  in  Banks  (Ref  15).  Since  the  results  are  presented  as  a fraction 
of  vacuum  fluence  in  mean-free-paths,  the  conversion  to  build-up  factors 
must  be  performed.  A comparison  to  the  build-up  factors  calculated  by  the 
moments  method  is  presented  in  Fig  10?.  It  can  be  seen  that  the  percent 
differences  for  60  keV,  which  is  the  only  energy  range  given,  are  from  1?% 
to  5256. 

Krumbein  (Ref  16)  using  SAGE  performed  calculations  of  x-ray  air  transport. 
SAGE  is  United  Nuclear's  spherical  Monte  Carlo  photon  code.  The  calculations 
were  carried  out  at  various  energies  from  12  keV  to  100  keV  and  out  to  10 
mean-free-paths.  Even  though  the  calculations  for  the  fluence  used  to 
calculate  the  build-up  factors  were  performed  using  the  density  of  air  at 
an  altitude  of  22,860  meters,  this  had  no  effect  on  the  build-up  factors 
since  the  energy  build-up  factor  Is  independent  of  altitude  when  It  is 
reported  as  a function  of  mean-free-path.  Comparison  to  Krumbein* s 
calculations  is  6hown  in  Figs  108,109  and  110.  P'or  20  keV,  the  percent 
differences  range  from  7* 5%  60  1 7%.  For  40  keV,  the  percent  differences 
range  from  18%  to  32%.  For  100  keV,  the  percent  differences  range  from  12% 
to  27%.  These  differences  can  be  accounted  for  from  the  different  cross 
sections  used  and  by  taking  into  account  the  error  produced  in  the  calculation 
by  the  moments  method. 

One  of  the  most  extensive  calculations  using  a Monte  Carlo  code  to 
solve  the  x-ray  air  transport  problem  in  air  was  perfromed  by  Shelton  and 
Keith  (Ref  17).  This  calculation  was  perfromed  by  the  HAT  code.  The 
results  are  presented  in  comparison  to  the  moments  method  calculation  in 
Figs  111-118.  For  20  keV,  the  percent  differences  are  from  -8%  to  -69%. 

For  40  keV,  the  percent  differences  are  from  -2%  to  -41  J%,  For  60  keV,  the 
percent  differences  are  from  -8%  to  -76%.  For  90  keV,  the  percent  differences 
are  from  3%  to  -?0%.  For  120  keV,  the  percent  differences  are  from  7%  to 
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-27*.  For  150  keV,  the  percent  differences  axe  from  5*  to  -9%.  For  200  keV, 
the  percent  differences  are  from  1%  to  -63*.  For  3^0  ke^,  the  percent  dif- 
ferences are  from  -2*  to  -96*.  The  differences  go  up  as  range  Increases  In 
all  cases  execpt  150  keV.  The  absence  of  smooth  curves  In  the  range  verses 
build-up  factor  semi-log  graphs  indicates  sizable  statistical  variation  in 
the  build-up  factors  calculated  by  this  Monte  Carlo  code.  This  can  account 
for  the  large  amount  of  disagreement  at  some  ranges  and  energies. 

Discrete  Ordinates.  Because  of  the  limited  amount  of  published  results 
of  discrete  ordinates  calculations,  only  one  comparison  of  the  moments  method 
calculation  to  discrete  ordinates  will  be  presented.  Dupree  (Ref  18)  used 
DTFXRAY  to  perform  a limited  calculation  at  100  keV.  Since  energy  fluence 
was  obtained  for  only  three  distances,  a table  of  results  along  with  the 
results  of  the  moments  method  and  the  percent  differences  are  presented  as 
Table  II.  The  results  are  in  great  disparity.  The  erratic  agreement  is  due 
to  either  the  limited  and  truncated  discrete  ordinates  calculation,  which  was 
a P-jSg  or  opposite  maximum  error  in  the  discrete  ordinates  and  moments  method. 

Bigelow.  Because  the  program  used  by  Bigelow  was  the  same  one  that  was 
used  for  this  study,  comparison  would  not  be  in  order.  But  using  the  same 
input,  the  exact  results  were  obtained,  showing  that  the  program  was  working 
as  programmed  and  any  error  calculations  performed  by  Bigelow  are  applicable 
to  this  study.  When  recent  cross  sections  for  energy  range  from  1 keV  to 
100  keV  were  used,  deviation  in  results  occurred  at  small  energies. 


Table  II 

Comparison  of  Energy  Build-up  Factors  for  100  keV 
Mean-free-paths  Discrete  Ordinates  Moments  Method  Tercent  Difference 
2 4.80  9.12  47* 

4 13.6  29-0  53* 

10  97.6  205.  52* 
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Purress  arri  S??Pe 

The  purpose  of  this  report  is  to  provide  a complete  set  of  monoencrgetic 
build-up  factors  which  are  used  to  simplify  the  x-ray  air  transport  problem. 

To  compute  these  build-up  factors,  a moments  method  program  was  used.  The 
program  PHOTDIS  was  chosen  and  run  on  a CDC  6600.  The  energies  considered 
were  from  14  keV  to  1000  keV.  The  ranges  considered  were  from  one  to  15 
mean- f ree- patns . 

Results  and  Discussion 

The  results  were  presented  in  a set  of  graphs  for  each  energy.  The 
intrinsic  accuracy  of  the  program  estimated  by  convergence  analysis  performed 
by  Bigelow  was  found  to  be  a maximum  of  20%.  Comparison  to  others  found  the 
results  of  the  moments  method  calculation  to  be  higher  than  some  Monte  Carlo 
calculations  and  lower  than  others.  For  10  mean-free-paths,  this  is  shown 
in  Fig  119*  To  facilitate  use  of  the  data  presented,  selected  energies  are 
presented  together  in  Figs  120  and  121 « 

Recommendatlor.3  and  Summary 

Because  of  the  limited  time  and  the  nature  of  the  program,  experimentation 
with  different  number  of  moments  and  the  number  of  polynomials  vised  in  the 
reconstruction  could  not  be  performed.  It  is  recommended  that  any  other 
extensive  study  that  computes  build-up  factors  using  the  moments  method  should 
vary  the  number  of  moments  and  reconstruction  polynomials. 
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APPENDIX  A 


Solution  _of  th9_Boltzratm-.Tro.TisPQrt_EQuatlQfi_by_  th^J'.oaentsJ^thQd 


Reduction  of  the  Polt"<r.arn_?ra.nsPQr.t  Equation 

From  Eq  (5),  the  Boltzmann  Transport  Equation  integrated  over  all  time 
for  spherical  coordinates  can  be  extracted! 

arhuSdii  + I3l!  + utF(r,E,|1)  « //  F(r,E,n')nS(u  ,E'-E)dE'd4* 

(ll) 

S(r,E,n) 

where 

H is  the  cosine  of  the  angle  0 in  spherical  coordinates. 

Since  the  virgin  fluence  can  be  calculates  without  resort  to  computers  by 
Eq  (l),  the  fluence  in  Eq  (ll)  is  the  scattered  fluence. 

To  expand  the  functions  in  Eq  (ll)»  a set  of  polynomials  is  needed. 
Legendre  polynomials  conform  to  the  requirements  for  expansion  polynomials. 
These  polynomials  have  the  following  properties! 

•+1 


/_!  W.O'**  - K 

fw  ■ £ 

n~u 


/*+ 1 

fn  mJ  f(x)Pn(x)dx 


(12) 

(13) 

(14) 


Therefore,  the  expansion  of  the  fluence  F,  the  Scatter  cross  section  ^ , and 
the  source  S in  Lengendra  polynomials  is 


F(r,E,^)  - £ ^Fn(r,E)Pn(h) 
n=Q 


4*(lV*UI>  - •£.  ^ °0M)P („  ) 

m*U 


(15) 


(16) 
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+ SA.C*.*)^)2^  - E S^r.DP^)^  (18) 

n=0  n=0 

♦ // (§  f'»(r,E')p>')2S1) (£  w' 

By  rearranging  the  scattering  inxegral(the  second  term  on  the  right  side), 
it  becomes i 


By  using  the  addition  theorem  on  the  double  Integral  above,  the  integral 
becomes 

CCfoW'  * 2 SiSkH  eo3[fc(«-0')]| 


Upon  performing  the  integrals 


tyrr 


K 


2n+l  -nin 

So  the  scattering  integral  in  Eq  (l8)  becomes 


fo  [(&  F"(r,E)2^)  (S  v> B)2^)  * P>)S™ 


dE' 


which  equals 

: % 

Putting  the  last  term  into  Eq  (18),  the  full  equation  becomes 


F E '„<*.*)  fS(^E)Pn(p)dE' 
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wrm 


+ £ U>n(r.E)P>)  X?  -ft  £ ^(r.E-^tE-E)  ^ Pn(„Wl9) 

n°U  0 n=u 

♦ f;^  ■.«*m»,w 

By  multiplying  each  term  of  Eq  (19)  "by  2rrP  Cu)  dtl  aild  integrating  each  term 
from  -1  to  +1,  the  following  equation  is  obtained! 

*F  (r,E)  r+l 


2n±i 
2 


■V-/t  *>>*>>*• 


00  . /*+l  o ®P  m/  4. 

+ lJ‘2Slpn(r'E)  J t Cini2>-~ — uV^r.E)  - S^E^r)  (20 ) 

+ ^E  Fx(r,E')^(E’-E)  dE' 

To  simplify  Eq  (20),  the  recursion  relations  for  Legendre  polynomials  must 
be  employed! 

. r 1 

(21) 


P»M  " 2^1  [(nthl)W^  + "W^] 


(i-u  ) 


2\  _JQ 


■ M [t»l)Pw,(u)  ♦ ■VjCu}]!  <22) 

By  substituting  Eq  (21 ) into  the  first  tktrm  of  Eq  (20)and  by  substituting 
Eq  (22)  into  the  second  term  of  Eq  (20)  and  perfroraing  the  integration,  the 


final  P equation  is  obtained! 

X 


. _i_  dFWr,E)  . tetiKx±al  Ini^L — 


2M  1+1 

2.6+1  dr 


2£+l 
,E 


(23) 


2f+l  dr 

’ ^2^  + " Sf(E)&(r)  + /;  ^(r,E’)hS(E’-E)dE' 

It  is  necessary  to  introduce  the  definition  of  the  moments  of  a function  fi 

M - f °°rnf  ( r ) 4irr^dr  (24) 

n •'0 
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where  M is  the  nth  moment  of  the  function  f.  To  obtain  the  moments  of  the 
n 

fluence  in  Eq  (23),  it  is  essential  to  multiply  each  term  in  Eq  (23)  by 
rn4prr2dr  and  integrate  from  zero  to  infinity! 


—2±1 
2 X+l 


Z*00  n aFP+i(r»E)  ? p n 3F.  (r,E)  2 

J0  r — "T,r  dr  ♦ J0  r — ,imTiz 

.<*U)Ut2)  £ r»  /“  r“  — — “4tfrEdr  (25) 

30  |2  /*oo  »oo 

ti*  f rnF,(r,E)4Trr2dr  - f I r"VsF  (r,E)4Trr2drd£  + / rnS.(E)&(r)4irr2dr 
Z J 0 * J Q •'O  ^ *'0  X 


Using  the  moments  definition,  Eq(24),  the  following  equation  is  extracted! 


_J±i 

2X*1 


“ - afV,(r,E) 


n 

r 

rt 


itl 


dr  ^ dr  * 2*41 


°°  3F.  (r,E) 

— £_J L d 


f r° 


4rrr  dr 


(26) 


+ ■CfrtllC&'tZ.l  ^ M + |A%!  ■ S, 

2X+1  £+l,n-l  2X+1  Jl-l,nrl  Zn  *Z 


<e>*X  v6dE 


where  M,  is  the  nth  moment  of  the  Xth  expansion  coefficient.  By  applying 
£n 

integration  by  parts  to  the  first  two  terms  of  Eq  (26),  they  become 

~ 2X+!(r*2)Me+l,n-l  " 2i+l^2^!i-l,n-l 
So  the  full  moments  equation  is 

[ 2Z+V  2'  2X+1  J x+l,nrl  [2£+l^2'  2X+1  J Vl,n-1 


+ n " SX(E}  * hit1' 


(27) 


'dE 


Rearranging  and  combining  terms,  the  recursion  equation  for  the  moments  is 

.. tj,  _ fE  ,.sM  (v)  jjtfcf+la  M + £±2Lz2rji  „ (28) 

* HXn  J0  4 MindE  SZ'h)  + 2X+1  KM,n-t  + 2*+l  £-l,n-l  '‘2b; 

The  integral  is  evaluated  by  numerical  quadrature.  With  Eq  (28),  all  the 
moments  can  be  calculated.  But  only  moments  where  X-n  is  even  or  ir n is  odd 
are  related.  Also  moments  where  Z-n  are  negative  have  no  physical  significance. 
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Bs&Qrurtrw?UQPL.gg  the  fiance 

Even  though  Eq  (2b)  Is  for  the  scattered  contribution  to  the  fluence, 
the  once  scattered  fluence  can  also  be  calculated  directly  using  (flef  Iil2) 
. STu(£-cosa)lE(E  -»E) 

rkr.E.u)  - J ®— 


8rr  r(slna)(sin  0) 


|«p^|lt(s0)gj{t  Ut(E)S^il)]|  (29) 


where 


F Is  the  once  scattered  fluence 


cos 


1 1 + (mt 


c‘/E)  - (n  cVE) 


>] 


0 for  x<0 

1 for  x>0 


u(x)  is  the  unit  step  function  - 

Eq  is  the  source  energy 

ST  is  the  total  source  strength 

m is  the  mass  of  an  electron 
o 

Therefore,  all  that  is  needed  to  find  the  total  fluence  is  the  multiscattered 
fluence.  This  fluence  can  be  calculated  from  the  multiscattered  moments  given 
by 


M™  - m - jr 
Zn  £n  r‘fn 


(30) 


where 


K®  is  the  multiscattered  moments 
£n 

is  the  once  scattered  moment  calculated  from  the  once  scattered 
^ fluence 

To  find  the  multiscattered  fluence  from  its  moments,  a moments 
reconstruction  must  be  performed.  The  first  step  is  to  define  the  expansion 
function  of  the  moments,  G,  which  is 


<-i)J 


(2.fr.t)inl 


(2.W22;Hr.-j;iJI  "2.2+2J 


(31) 


Using  the  available  moments,  the  expansion  function's  values  can  be  used 
to  calculate  the  expansion  coefficients,  F*!|(r,E)  for  the  multiscattered  fluence  1 


148 


F%:.£)  - f]  a,  V 


where 

“wl-zfe  (a.+  W+l-yJwJ+j^  (33) 

and  ® 1.  It  should  be  noted  that  in  Eq  (31)*  only  the  moments  where 
X-n  are  positive  and  even  are  needed  to  calculate  the  expansion  function. 

The  nultiscattered  fluence  is 

Fj(r,E,n)  = Fj(r,E)Px(d)  W 

So  the  total  fluence  is 

F(r,E,n)  - F°(r,E,p.)  F1(rtEp|i)  + Fm(r,E,u)  (35) 

where  F°(r,E,d)  is  the  virgin  fluence  calculated  from  Eq  (l),  F*(r,E,d) 
calculated  from  Eq  (29)  and  Fn(r,E*n)  calculated  from  Eq  (3^)«  The  build-up 
factor  is  given  by  th<  following  equation i 

lU  (36) 

F°(r,E*n) 


1 


APPENDIX  B 


Sample  Input  and  -Program.  .LlaUng 

Sample  Input 

Th®  sample  Input  listed  on  paces  151"15^  is  the  input  used  in  this 
study.  This  listing  is  the  data  input  into  unit  5.  Unit  4 inputs  the 
number  of  energies  to  run,  the  energies  and  the  ranges  in  mearrfree-paths. 

Prorratn  Listing 

The  program  listed  on  paces  1 5^"  177  is  the  version  of  PHCTDIS  used  in 
this  study.  This  program  is  exactly  the  same  as  the  program  used  by  Bigelow 
(Ref  l)  execpt  for  format  of  input  and  output.  Output  is  on  units  2,3*  and 
6.  Unit  2 is  output  for  input  of  a plotting  program  used  to  generate  the 
plots  on  pages  16-120.  Units  3 and  6 are  printed  output.  Much  intermediate 
output  is  eliminated  for  this  study,  since  fluence,  moments  and  convergence 
analysis  are  not  desired.  This  program,  run  on  a CDC  6600  takes  48  central 
processor  seconds  and  260  input-output  seconds  per  energy  using  five  space 
points. 


150 


NIT  9 o'* E N 

7 

. 2. 

c 

O.iOODC'O^E- 

33 

0 . 

■’OOOOOPE 

01 

0. 54-'  0 0 0 0" 

00 

0. 

5 3 'i  9 0 0 0 F 

01 

0.1 3700m" 

01 

0 . 

3 3 P i 0 0 0 E 

01 

0.19^0  0 0 3E. 

Cl 

0 . 

n^onooE 

01 

0.?*jF  3 00  0E 

01 

0. 

l n n o o n c 

Cl 

0.131009 0E 

01 

0 . 

1 o 3 o o 0 E 

01 

0.3  90  0 0 DOE 

01 

0. 

1 ° 3 3 0 0 OE 

01 

0.4  V:  OOOQE 

01 

0. 

iooo  no oe 

Cl 

0. 531000 QE 

01 

0. 

'■•ooo  oooc 

OP 

0 . :>  97  0 9 9 JE 

01 

0. 

= 000  ? 0 0 E 

00 

0 . f-  5 1 0 0 9 0 E 

01 

0. 

3 0 : 1 C 0 0 E 

oc 

0./300000E 

01 

n. 

5 TOO  DOPE 

oc 

0 . 12.;  ICO  9E 

(V’ 

0. 

?7  v 3 0 0 E 

or 

o.eooocooE 

02 

0. 

io  no  n ooe- 

01 

0.7 000 00 OE 

02 

-0. 

30r  0 o ooe 

00 

0 . 1 0 0 0 0 0 0 E 

03 

-0. 

OOP p 300E 

cc 

0 . 1 0 0 0 0 0 0 E 

3? 

-0. 

0 9 9 n 9 0 05 

on 

1.0  00 

7.  9 53  OO  r 0 4 

2^0  0E- 

Cl 

1.500 

2.enr.oor  ^ 

F. 15700E- 

01 

2.000 

l.ll  -’00"  0 4 

t .1 239CE- 

01 

3.000 

3,4  4 3 0 0"  0 3 

1 . 7 8 1 0 0 E 

00 

4.000 

1.4  70  OOF  04 

1.6  54  OOF 

00 

G.COO 

4.n0’99"  02 

2.  *'  8 700 E 

30 

6.000 

1.5-04  9 0 c 02 

?.*  EfoOOE 

00 

10.000 

’.7  31  .DDE  01 

3 , C 9 5 0 0 E 

90 

15.000 

2.  <U  “ 3 0r  0 1 

?,.i  F OC E 

09 

20.000 

r.«-oc'oor  oo 

3.1'  51  OOF 

00 

30.000 

2 . 0 •"  r<  0 9 r 0 0 

3.7  0200E 

00 

40.000 

7.517  COE- 01 

3.7 1 1C0E 

00 

60.000 

2.074305-01 

3.'  9900  E 

09 

30.000 

P .10’00F-02 

3.:  4900E 

00 

100. OOC 

3.? 77 00 2- 02 

7.1  0 2 COE 

90 

150.000 

i . ore  oor- 02 

3.  f c 3 0 0 E 

00 

200.000 

4.40OOOr-O3 

2 ,f  3 3 00  E 

DC 

300. OOC 

1.2  6*0^-0  3 

2.  tfcf  03E 

90 

400.000 

5 .5  19  OOF-  P*i 

2. 21300E 

9 0 

500.000 

2.5 17 30 r- 04 

2. C2200E 

00 

600.000 

1,75"  0 3 c - 0 

l.S 7 100E 

00 

600.030 

B . n-,3  0 0r-  05 

1 . f 4 4 0 C E 

00 

1000. 000 

l . 3 2 7 00  E-  35 

1 ,i  7 8 00  E 

00 

0 . 1 4 9 3 6 E 

04 

0.?"009r-94 

G.12G0FE 

01 

0 . 1 5 0 0 0 E 

04 

n.250nnr-n4 

0.120  OF  E 

01 

0.20000c 

34 

0.1r0O0r-O4 

0 . 1 0 23!'.  F. 

01 

0.  339  0 OP 

04 

0 . 3 P 0 0 0 F - C 5 

G.(  0 ^ 6 3 E 

00 

o.noooflp 

04 

0,7r'000r-0  5 

O.i"  7 064" 

0 0 

0.50000E 

34 

0.2  30  00r-05 

n.‘  7675  E 

40 

P.bOOOOE 

O'. 

0 .2  ‘0  mr-03 

0.f  1121 E 

no 

0.80900E 

Oc. 

0,150 03c-O6 

0 • 1 1 7 6 7 E 

00 

0. 100  o<;p 

05 

o.i or  oor- op 

0.  3f 5 32 E 

00 

0.1 5 OTOE 

0" 

n.90000r-C3 

u . ? 6 ? 2 d F 

30 

0.20000E 

0« 

0.5f000r“0b 

0. 2 0 9 9 j E 

09 

0.300001 

05 

0 . 1 0 0 00r‘-^0 

L.  15211  E 

OC 

151 


0.4 0330^ 

Or 

C.1  3p>nir-^0  0 

. 1 2037E 

00 

0.5  00  3 0'" 

Or 

0.lonoor~70  o 

. 1 00 J2E 

30 

0.600006 

0"> 

0. 10030c- 70  0 

. f 6?3i  E- 

01 

O . ft  0 3 0 0 r 

04 

0 . 103  OOP- ^0  L 

. i.7  19?6- 

0 1 

0.10330': 

05 

0.1000 3 r-30  0 

. f 64  0;  E- 

01 

xvo:‘i 

8 . 2 . 

6 .209**7 

(.  OE 

0.  ICO  3 0 3 0 

r- 

0.  3 30°  0 OOE 

01 

0,r>0»-»3  03  0t 

00 

0. 64 ooo ore 

01 

0.1373030 

c 

01 

o. 7o  on  o orr 

01 

0.1930000 

ZT 

01 

0 . ?•-  0 ')  0 0 0 r 

0 1 

0.3530030 

r 

01 

0 . 20  0 3 3 0 OE 

01 

0 . 3 3 1 0 r .3  0 

c 

01 

0. 1603 n OOF 

Cl 

0.33^0000 

TT 

01 

3 . 1 - n r*  0 n C- 

01 

0 . ft . 3 0 3 3 

T 

01 

n . 1 3 o •>  o c c r 

01 

0 . •>  3 1 0 0 J 0 

01 

o.  rn  poof 

01 

0.9970  C OO 

r 

01 

0.  110”  ,,00E 

01 

0.6C'-3033 

c 

01 

0. 60030005 

or 

0.7 30 0330 

01 

0.8303000E 

or 

0 . 1 2o  1 0 0 0 

T 

OP 

0.  r-  3 03  c 0pr> 

01 

o.pooonjn 

C 

02 

- 0.  03  0 ) 0 0 06 

or 

o. fooocoGE 

C? 

-0.  3 30  3 hope 

or 

0.1003030 

r 

33 

-0.  000>3  300F 

00 

0.100  3 OTOE 

3? 

- 0 . 03  33  n 0 CE 

or 

1.0  00 

1 . 2 2 1 0 3 r 3 5 J 

.1  7 300E- 

01 

1.-300 

4.157 00 r 04  3 

170CE- 

c 1 

2.000 

1.373006  04  6 

.7  C*-Orc- 

0 1 

3.003 

3 .705  00  6 os  1 

. 220  0CE 

00 

4.000 

2.  *4  51  0 36  0 7 1 

. 7 7 n CO  E 

00 

5.000 

1.240 03r  33  2 

. 1 7 1 0 0 E 

30 

6.0  00 

7. 045  0 IE  3 2 2 

. r 40  00  6 

30 

6.003 

2 . S 3 i 03r  3 2 ? 

. C6300E 

CO 

10.030 

1.400036  0?  3 

. 754 00 E 

•3  0 

15.000 

3 . 7 3 ^ 0 3F  31  3 

,r4403E 

30 

20  . 00  0 

1.4  31  00  6 0 1 4 

. rsr  ooe 

0 0 

30.000 

3.921 n0  F 00  i 

. 27  ft OOE 

00 

40.000 

1 ."25006  00  i 

.709 OOE 

00 

50.000 

7.  ?r  0 3 36-  0 1 4 

. 27  1 03 E 

00 

60.000 

4 .054  00c-  01 

. 20  3C0E 

00 

50.000 

1 .6.96  Cnr_  01  L 

. 04  20PF 

00 

100.003 

7.7?? 306_n2  3 

76  OOE 

00 

l*;  o.ooo 

2.173036-32  3 

.l  1 9 OOE 

00 

200.000 

3.8  07  006-  3 3 7 

. 2 3 6. 0 0 E 

00 

303.000 

2. ”6 0036-07  2 

. 6 2 0 0 C E 

30 

4 o o . o n o 

1.10'»39c-33  2 

. r 2E0C6 

jO 

500,000 

5.?"*?0Or-04  2 

. 7 1 0 OGE 

00 

600.3  0 0 

3.522  006-  0*.  2 

. 1 ?eocE 

00 

6 0 3 . 0 0 (1 

1.5  15  33  6-  0 4 1 

.1  7 8 0rE 

00 

1000.300 

9.0?5  30  r-3  5 1 

.1.  6° OOE 

00 

0. 150  0 GE 

04 

0.450006-34  0 

. 1 37  21 E 

01 

0 , 20C  3 0r 

Oi 

0 , 2 0 0 9 3 r - 3 4 G 

. 1 1701r 

01 

0. 70000E 

04 

0.1 2" 306-0 4 0 

.c  i«5?£ 

00 

0 , 4 0 0 0 0 E 

04 

0 .P"  0 03E- p5  C 

.7604 4 E 

30 

0.50000F 

04 

0.5 20 OO6-0 E c 

. I 6 1 1. 3 E 

00 

152 


0.50  3 0 1 ' C O'.  0.:- 30  00r- n">  G.'6‘.?i,5  U 0 

0.530009  O'-.  0 ♦ ' 7 7 3 , F 00 

o.iooooe  oi-  o .2n  noc-'i5  i.t  geu7 e oo 
0.150109  05  0.13nmr-u'.:  0 . 2 9 ° 7 2 E 00 

O.FOOOO17  0"  O.no  115-00  G.?39"JhE  00 

0.300  0 09  oK  o.innm^-^Q  o . 3 7 ?e-  =:  00 

0. < 00  0 09  03  0. 10  0 03— ’0  0 . 1 37 1-  ? 9 00 

0.500005  0?  0.1011',r-’0  t.  11031'*  00 

o.50oanr  o“  o.ionoE-n  G.50O95E-O1 
0.600005  OS  0.100005-70  0.7  6791E-01 


ARGOM 


13 . 1 . C . 93  t OOOOE-O? 


0.10000105- 

3C 

0.  13  no noo 

99 

02 

0 ,9'j:  0 0 0 09 

00 

n.  r goo<7 

0 2 

0.13^00009 

01 

o . 1 0 7 ' 0 0 

r 

02 

0 . 1 j 1 0 0 0 0 - 

01 

0 . 1 2 0 0 1 0 c 

c 

Cl 

0 . 2b:i  0 0 0 09 

01 

o.'-onnooE 

01 

0 . 331  0 0 0 19 

01 

0 .33  ncoo 

c 

01 

0.33’,0  CO  09 

01 

0. “30 T 10  0 

r 

Cl 

0 0 0 0 03 

01 

0 . n 0 j z 0 n 

c 

01 

0.9310000- 

01 

0 . 3'n*  o 0 09 

01 

0.9 9 7 00009 

01 

P.’lOnpGO 

c 

01 

o.p.^'+nooo9 

01 

9 . moo 

r 

01 

0.73000009 

01 

0.  ’201100 

r 

01 

0.1?Mn09 

0’ 

C.  nn  100 

c 

01 

0.  ’bn  00  r 9 

0’ 

0. 3001000 

c 

on  • 

0.70000009 

02 

- 0 . nnioon 

r 

00 

0.10000009 

03 

-0.  n on 0 0 

C 

0 0 

0 . 1 00  0 0 0 09 

3? 

-0.  0 0 0 n i 0 0 

9T 

On 

1.000 

? 

. l.noo-  05 

.(  17  005- 

f 1 

i.r-0  0 

7 

. 0 3 0 0 r 0 if 

O 

. 33*  005- 

0 1 

2.  000 

3 

.40 "in  nn 

j 

.'55005 

0 0 

1 3.000 

1 

. 1 1 0 0 0 r 04 

c 

.<  5 3 00  5 

00 

3.  ’01 

9 

.297009  03 

2 

.f-  ??O0E 

nC 

3.20? 

0 

.165  309:  04 

2 

.4  27005 

30 

'4.000 

* 

5 

.223105  Oh 

3 

• 322005 

OC 

5.000 

2 

, 9 09  n0r  04 

1, 

. C 39  00  5 

00 

8.000 

7 

.9  7 2 0 0 F 0 3 

r 

• ! 11 0CC 

00 

10.010 

4 

.2  03  00-  0 7 

6 

. 1 62005 

00 

15.000 

1 

.2“ 0305  03 

■1 

1 

.35  200  5 

00 

20.000 

5 

.7  mi r 0 2 

b 

. 0 7 0 00P- 

00 

30.000 

1 

. 5 ? 2 00r  C 2 

0 

.7  4 8 0 C 9 

00 

*-0.00  0 

5 

. 2 0 r 0 0 9 n 1 

s 

.5 60  ODE 

OG 

J 50.000 

3 

. O-HO lr  "1 

9 

. r 40005 

0 0 

eo.opo 

1 

. 74100r  01 

p 

.9  99005 

10 

60.000 

7 

,07700r  00 

t 

.76  0 009 

00 

103.010 

3 

."•23oor  n' 

r 

* • 0 0 0 G 9 

00 

15  0.000 

1 

,012006  00 

7 

. l-  0 3 0 1 c 

00 

200.000 

'4 

.2  39  .115-  0 1 

7 

.206009 

00 

300.010 

1 

.26301 r- 01 

r 

U 

.311009 

00 

A 0 0 . 0 1 0 

5 

.644 10  r-  1 2 

U 

.C  71006 

0 0 

900.000 

3 

,04’00r-C2 

4 

. 1 6 b OG  5 

00 

000.010 

1 

.or  n nr.  <ig 

U 

.*  57  009 

00 

300.030 

p 

v 

.500016-03 

L 

. 223015 

10 

1000. 000 

r 

. 040 005- 03 

3 

.797  009 

00 

kt: 


F • — - — ----- 


0 . 1 - 0 0 c E 

O'* 

0.?7‘,0  0r-  1? 

0.  703!  ?? 

01 

D . ? 0 0 D C £ 

0-t 

0.1  "7  '’0  — 0 2 

0.2U217E 

0 1 

0 . 3 1)  0 0 II E 

ov 

c.9f-' oo—o.i 

0 . 7 C 66 3 E 

01 

0 . '•  0 0 l)  0 c 

04 

0 . P- ? 9 1 c - 0 3 

1.17  2 1 • E 

01 

0 . 5>  0 0 0 P E 

O'. 

0 . 1.  V->  00 r - 0 3 

0 . 1 466?17 

01 

0 . r>  0 o a 0 r 

3', 

0 . 3 ^ 0 0 - - ''3 

G . 1 314;F 

01 

0.3003CE 

04 

0.2. 01007- 03 

0.  1 0 7 <*  0 E 

0 1 

0.3  O'}  OPE 

00 

0 . ?fl 0 0r-0  3 

0.1  0 7 <-  C E 

01 

O.OOOOOE 

C-. 

o.?r,inr-o3 

G.  1 07  /■  Or 

01 

0.303 3 CE 

00 

0 . 0 31  0 0E-  ■"» 

0.1  07-.0E 

p 1 

0 . 5 0 0 0 0 r 

C4 

o.?pi  n--'13 

0 . 1 07  4 OE 

01 

O.oOQOOE 

04 

0 . 2 r 1 0 0 r - ' 7 

0.1  0740  E 

PI 

0 . * 0 0 1 P E 

O'- 

P . ?f 1 on-- 0 3 

0 . 1 0 7 ^ 0 r 

01 

0.300 3 CE 

O 

0-t 

0 . ? S 1 00  E - 0 3 

0 . 1 0 7 4 C E 

01 

C 

0 . 0 00  0 PE 

00 

o.m.’ie 

20 

0. 200  3?E 

3 ^0  1 
ti  21  1 

03 

b 6 13 

o.^cnnnooF  24 


0.0033  41 

91  nr- 

r 

61  111 

331 

1'  21  31 

4 1 

91  371 

1.00  0.99 

0.90 

0.=: 

0.94 

0.12 

0.90 

0.06 

0.56  o.n 

0.6?  0.60 

...  0 . y 0 0.60 

0 . 40 

0.0  0 

-0.20 

-0.'  0 

-0.60 

-0.50 

- 1.00 

15*» 


1 

1 


oooo  ooo  o o ooooooooooooooooooo 


a 


PROGRAM  PHOTO  IS r NPUT , OUTPUT , TAPED, T A PE6 , T ft  PEC , T APE3 , T ftp E4 , T ftp E3 , 
1TAPE2) 

SUBROUTINES  USED  BY  THIS  PROGRAM  SS$$$$SJS$ 

SIG CAL DULft Tr S TOTAL  MACROSCOPIC  CROSS  ScCTION 

ft  KERN CALCULrTES  DIFFERENTIAL  SCATTERING  CROSS  SECTION 

(SCATTERING  KERNFL ) 

POLKER CALCULATES  PPODUCT  OF  LEGENDRE  POLYNCMIA,  AND  SCAT- 

TE  °I N"  KTRNFl 

ONCE  r V A L U A ~ R S PART  OF  THE  ONE-SCATT  rR  FORMULA,  rD  (?1) 

KXPAN IJFrS  P"runSIO(  RELATION  TO  CALOULATF  COErF 1 0 I E N TS  rOR 

THE  EXPANSION  POLYNOMIALS  ( EO  (DA)  > 

BAITA CALCULATES  EXPANSION  COrFrICT~NT  OF  EO  (27) 

FIRST CALCULATES  ANGULAR  DISTRIBUTION  OF  ONCE-SC 5 T T ERFP  ?HO 

TONS,  "LSO  ALL-ANGLE  FIRST  SCATTER,  AND  rI  7ST  SCAT T FR 
ENERGY  DENSITY 

SPACE I NT  ERL  T MK AG F OF  THE  MOMENTS 

INC CONTROL0  INDEXING  OF  ARRAY 

TORIAL CALCULATES  FACTORIALS 

COMMON  PHOT  (3  , AO)  , ENERG  ( 7,4  0)  , INCO ( 7, 4 0)  , °T A (7,  17) , CAPK(  3, 1 7)  , ON ( 3 
1)  , DN7  ( 3)  ,NXS,NrL»T0LF,TCLK,PIE,R7ERD,LMftX,  XMU'-'U  , X MUL  A M , T EMP,  YlY,  DM 


TORIAL CALCULATES  FACTORIALS 

COMMON  PHOT  (3 ,4ft)  , ENERG (3,40) , INCO ( 7, 4 0)  , FT  A <3,  17) , CftPK( 3, 17)  , ON(3 
1)  , DN7(  3'  , NX  S , NrL  , T OL  F , T C L K , ? I E » R7  ER  0 » L M A X , XMU'-'U  , X MUL  A M , T E HP,  YlY,  DM 
?U , SIGMAG , SI  NT  ML , ~ I GO , ONE , XLNCK  > L , N, L P , N p , XL  AMG , BET  A , 5 ZERO , 

3NTRANS  ,NHU,XL  AMO,  GAPB,NYP,NCUT  , FACT  0 R , BSC.  AT  ( 10,  19)  , 

5 SIG  MAT  (331)  , L AMBB  A ( 33  1 ) , ONESCT  ( D , 3 7 , ? 1)  , AL  ANG1  ( 5,33)  , ALFNG  (5)  , XMU( 

6 21)  , Y ( D)  , p M !J  ! 2, 21 ) ,VIRG1 N (D) , ALENGM (5 ) , A ( 10 , 19) 

DIMENSION  AnR  AY ( 1 0 , 19 , ?E)  , DEL T A D( 10 , 1 B)  , OMEP  (E , 33 , 21 ) , SU MN ( 1 0 , D , 33 
1)  ,pC(  1 0)  ,P9(  1 O)  , p L N K ( 1 0 ) , PL  NRM  1(10)  , PLNKM2  ( 1C)  , PL  N1'  *1  7 ( 10)  ,7(3)  , N AT 
30 MS (3)  ,WORD(7,?),NC(9),D(9),BONE(10,1B),INTEG(10,1B),B3(1U,19>, BSM 
41  (10, 19), PS  M2  (10,  19) , BSM 3 (10, 19)  , BON F p ( 1 0 , 19 ) , NSUMND ( 

510)  , ICOM(  16) 

REAL  L AMBDA , T NCO, N ATOMS , INTEG 
INTEGER  G , G 1 , G2 

DATA  CLITE/2.997925E10/,  FORPIF/ 12.55537/ 

RECURSION  RELATION  FOR  LEGENDRE  POLYNOMIALS 

ALEGRC  (X  , ARC,  POLY  1 , POLY  0 ) = ( ( 2 . * X + 1.  ) ' ARG*  POLY1- X*PO'_  Y 0)  / (X  + l.  ) 


»*»*♦♦»»»  <•**»»**>  K 

READ (4 ,1007) JOBF 
W R I T E ( 3,  201  7) 

WRITE (6, 2013) 

DO  999  JOP=  1,  JDBC 
F ACT  OR  = 5 1 1 .006 
DO  996  NSTU  F=  1,  1" 


nsumnq  (NSTiir)  =i 
AtPHA=  1. 

NPOINT  =1 

input  cross  srcrrnN$  fof  nit rogfn , o * y gen  and  argon-  eta  and 

n AP  K ARE  PARAMETERS  USED  TO  DETERMINE  HE  0 OUNO -ELEC  T RON  CCRRrCTIO 
TO  THE  <LEIII-NIHTNA  DIFFERENTIAL  CROSS  STCTION  FORM JL A 

00  150  1=1*3 

REAO(E>  ,10  10)  ( WOR''  (I,J),J=1,2),7(I),NATDMS(I),0N(I) 

READ (5 ,1  000)  (FT A ( T ,K)  , C A P K ( I , K ) , K = 1 , 1 7 ) 

READ  (5 ,1000)  ' CNF?G  ( I , l< ) , PHOT  (I  , K)  , I NC 0 ( I , K ) , K=  1 , 4 0) 

CONTINUE 

I PASS=  2 WAVELENGTH  MESH  PARAMETERS  TO  PE  INPUT 

= 1 OTHER  NT  S r 
REA  D (5 ,1007 ) ir>ASG 
ALTITUDE  AND  PARTICLE  DENSITY 
REAO(5,100F)  ALT 
H R I T E ( 3 , 1 100  ALT 
R F A 0 ( 5 ,1.003)  OFMSTY 

ENERGY  A GOVE  WHICH  NO  CORRECTION  IS  TO  G E HADE  TO  THE  <L  FIN- 
NISH I N A SCA TERTMG  LAW 
READ (5 ,1 003) c NX  N 
XLNCR- FACTO  R/ENRN  * . 00  01 

JND=2  FOR  PHASE  II  (RECONSTRUCTION)  ONLY 
=1  OTHERWISE 
RFAD(E»  ,1001  ) NFL  , XS  , JND 
NYP-l.G  ...  NUMr>rR  OF  SPACE  POINTS 

LEND=1  ...ANGULAR  PHOTON  DENSITY  IS  TO  PE  CALCULATED 
= 2 CALCULATE  THE  ALl.-ANGLF  DENSITY  ONLY 
R F A D ( 3 ,1001  ) NY°,  *1  MU  , L END 

NSTOP/  2=  1)1)  MG  - R ?r  MOMENTS  TO  PE  USTD  IN  I TT  ALL  Y IN  RCC  ON  STRUC  T T O') 
OF  THE  ALL-ANGLE  DENSITY  . . . N A D D / 2 = MJMRFR  OF  MOMENTS  TO  GF  Amro 
TO  THE  RECONSTRUCTION  EACH  TIM  F THROUGH  U>'TIL  N M A X 7 2 IS  REACHED... 
REA  0 (5 ,1001 ) NSTO~, NADD, UMAX 

E7FRO  = SOU^CP  F'JPRGY  *****  Q 7.  FRO  = SOURCE  STRENGTH 
R FA  0 (?  , 1000)  OFrpp 
R F A P ( A ,1100)  E7cRO 
IF (EOF  (4)  . NE.  0. ) GO  TO  DPO 

**J«*#J**V.<M*******»f**IMM  ft* 

$$$$$  UNITS  ...  'N  GENERAL,  CGS  UNI TS . . . FNERGIES  IN  KEV 

1 NO"  1 
LPJ  = 1 

P TE  = 3 . 14  139 
R7ER0=7.939BE-?6 
TOLE=.  00  0 1 
TOLK=.  00  0 1 
LMAX-NMAX/2 
L M A XP  = L‘l  AX  * 1 
NMAXP’-NMAXM 


oooooo  ooo  oooooooooooo 


NST0PP  = NST0P«- 1 
00  72  1=1, URL 

ON ( I ) = DE NST  v*  ON ( T ) *NAT0WS(I) 
0N7 (I) =0N (I) * 7(1) 

00  72  K=l,17 

CAP  = 1 , -CAPK  ( T ,K) 77  (I) 

IF(CAP)1  31  , 131,  132 

131  C A P K ( I ,0  = 0. IE- 30 
GO  TO  72 

132  CAPK(I ,<) =CAD 
72  continue: 

XLAM0=FACT0R/F7FP0 
G = 1 


DETERMINE  THE  WAVELENGTH  MESH  TO  BE  USED  IN  CALCULATION  OF  MOMENTS 

NTRANS  STE?  SIZE  SHIFTS  FROM  .02  TO  .04 

NOL  . . . SHIFT  FRO  H .04  TO  .00 

N0LP24  ...  ALL  MOMENTS  FEOUIRED  FOR  INTEGRATION  CAN  3E  STORED  IN 
ARRAY  (MOMENTS  F0D  20  VALUES  OF  THE  WAVELFNGTH) 

NCTOT  ...  1+NUMRer  0F  INTEGRATION  STEP  SI7E  CHANGES 
NC'.JT  ...  NUMBER  OF  MESH  POINT  AT  THE  FIRST  SCATTER  CJT  OFF 
MOC'JT  ...  NUMBER  OF  MESH  POINT  WHERE  SECOND  SCATTER  CUTS  OFF 
NBIGGP  ...  NUMBF5  OF  POINTS  IN  THE  WAVELENGTH  MESH 
NC  ...  LOCATIONS  OF  STEP  SIZE  CHANGES 

GO  TOT 40 3, 7 91) , IDASS 
403  0 ( 1 ) = . 02 

NC  ( 1)  = 1 
• NC  ( 2 ) = 21 
NTRANS  =21 
NC(3)  =71 
NC(4)  = 331 
N C T 0 T = 4 
NDL  =7 1 
NDLP24=95 
NCUT=6 1 
NOCUT=  91 
NBIGGP=331 
GO  TO  404 

ft*************#*-'  «**•**¥<*■»#»¥¥*****,  *¥***''*¥*¥*tH-**Mif*t  **¥¥***»*¥ 

791  READ (5,1012)9(1), NTRANS , NOL , NDL P24 , NCTOT ,NCUT, NOCUT, NBIGGP 
READ (5, 10 11 ) (NC(T ) ,1=1, NCTOT) 


SET  WAVELENGTH  VALUES  FCR  THE  MESH  AND  CALCULATE  CORRESPONDING 
TOTAL  CROSS  SECTIONS  (CM**-!) 


4 04 


10=1 

HCC-NC  (?) 
DFLT  A=  0( l ) 


157 


1 

i 


NRIGG=NBIGGP-i 
NGMAX=N5IGG/10 
LAM9DA (1) =XL AMO 
NGG  = 1 

71  ENERG Y =“ ACT  OR /L  A^  nDA (NGG) 

CALL  SIG(EMEPGY,GTG»-4AT(NGG)) 

I F ( NGG -N BIGGP ) 7 0 , 6 9 ,69 
70  NGG=NGG*1 

IF(NGG-NCC) 406,406,405 

405  f)ELTA=PFLTA+nEL7S 
TD=IO*- 1 

0 (ID)  = DELTA 
NCC  = NC  (ID+1 ) 

406  LAMBDA  (NGG) =L AIR ha (NGG-1 ) +0ELTA 
GO  TO  71 

69  CONTINUE 

CALL  PO.KFR(°C,XLAHO,LAMBDA(NCUT)  ) 

SIGO-SIGMAT (1) 

SIGOSQ=3IGO''SIGO/FORPIE 

SIGOIN  = . 01/SI  GO 

WRITE  ( 3,  1 1 0 2)  E7E?O,SIG0JN 

k 

9 

I F ( JND  .EO.  ?.)  GO  TO  126 

k 

9 

WRITE (6)  NRIGGDtLMAXP,NMAXP 
WRITE(9)  NBIGGD,LMAXP,NMAXP 

? $ S $ 5 $$$?$ 

BEGIN  PIUS’"  I 
S$$$?  $?$?$ 

NGG  = 1 
MCTR=2 
NOUT=l 1 

ASSIGN  734  TO  TGO 
ASSIGN  783  TO  ITRAP 
ASSIGN  259  TO  IFOR 
GO  TO  17 

RETURN  3 0 1 N T FOR  REGINNING  OF  CALCULATION  OF  MOMENTS  FOR  THE  NEXT 
MESH  POINT 

361  J=1 

LAST=1 

DO  35?  L P= 1 » L MA X° 

DO  352  NP=1,NMAXP 
352  INTEG(I.P,NP)  =0. 

XLG=L A M90A ( MCTR) 

NGG=MCTR 

GO  TO  ISO,  (761,784) 

784  JUMP=1 

ASSIGN  267  TO  NSK 
NOISCO=l 
L S E G = 0 


* 


158 


ooo  oooo  one. 


NSKIP= 1 


PIK'D  K MET  HP?  01 P"  ONT INUI T Y AT  FIRST  SCATTER  CUT  OFF  IS  WITHIN  THE 
RANGE  OF  INTEGRATION  AND  LOCATE  THE  LOWE'3  LIMIT  OF  THE  INTEGRAL 

I F ( MCT R .GE.  NOCJT)  GO  TO  300 

301  NGS= 1 

IF(MCTR-NCUT)  30  5,  3 05,303 
300  JUMP=2 

XLAMG2=L AMRna (MCTR) -2. 

DO  3 04  NGS-Nf'UT,MCTR 
TEST  = XLAHr,?-L  AM??  A (NGS) 

IF(  AOS (TEST)  -l.E-  05)  3 05,  3 05,  610 
610  IF(TEST) 305 ,305, 304 

304  CONTINUE 
303  NDISCO=2 

XL A MG2  =L  AM9QA (MC'R) -2. 

NGS=1 

IF(XLAM3  2-XLAMO)  7 05, 3 05, 313 
313  DO  307  NGS=2,NCU~ 

TEST  = XL6MG?-L AH50A  (NGS) 

I F ( A RS (TEST) -l.E- C5)  3 05, 3 05,  312 
312  IF(TEST)  305 ,3  05,7  07 
307  CONTINUE 

305  L IML  OW  =NGS 
GO  TO  302 

306  LIHLOW=NGS-l 

302  NBKSP=N3G-LIMLOW 

POSITION  T AaE  9 T0  PEAD  IN  MOMENTS  ...  LOOKING  FOR  THE  MOMFNTS 
CORRESPONDING  TO  THE  LOWER  LIMIT  OE  THE  INTEGRAL 

I F ( NRK SP  .GT.  LIMLOW)  GO  TO  750 

751  DO  752  I=1,NRKSP 

752  BACKSPACE  9 
GO  TO  753 

750  REWIND  9 

READ(9)NRIGGP,LMAXP,NMAyP 

753  READ (9 ) NSS, ( ( RS ( L P , Np ) ,LP=1,LMAXP) ,ND=1,NMAXP) 

I E ( NSS  .NE.  LTMLOW)  GO  TO  753 

CALL  PO'.KER(dO,  VL  AMO,XLP-) 

I F ( LIMLOW  .NE.  GO  TO  354 

DO  355  L P = 1 , L MA  X ° 

355  PLNK(LP)  =PO(LP) 

GO  TO  310 

354  CALL  POLKER(PLNK, LAMBDA (LIMLOW) ,XLG) 

FIND  WHERE  DELTA  LAMBDA  CHANGE  IS  THAT  IS  CLOSEST  TD  MCT  R 

310  DO  502  ID=1 , NOTOT 

IF(MCTR-MC  (ID))  717,317,5  02 
502  CONTINUE 


c o o o 


siinoiurar  7"r  rN’EGPAL  SO  THAT  SEGHrUTs  A M E FORMED  W H I ' H HAVE 
NO  CHANGE  IN  STEP  SI7F  WITHIN  THEM  A NO  NO  0 I SCO  NT  IN  JI  T Y 
LNOSEG  IS  THE  MESH  POINT  AT  THE  END  D E SUCH  A SEGM^T 

317  NCC  = NC  (I D - 1) 

IF(NCC-l)  73E,  736,  33.9 
336  L NDSEG  = MCTR 
GO  TO  555 

319  IE{  LAMBDA  (NGG  -1 ) -LAMBDA  ( NCC ) ) 264,261 , 262 

262  IE (LAMBDA (NGG -2) -L  AM3DA ( NCC) ) 264,  26  3,  26  4 
261  NSK I P=  3 

LSEG  = NCD-4 
ASSIGN  ?66  to  NS< 

GO  TO  572 

263  NSKIP=2 

L SEG-N  CD  - 2 
ASSIGN  266  TO  NSK 
GO  TO  57 2 

264  NSKIP-1 

ASSIGN  267  TO  NS< 

572  NPOINT  =1 
J = 1 

DO  562  ID=l,NCTOT 

IE(NC(  ID)  -LIMLDW3  5 6 2,56  2,563 

562  CONTINUE 

563  NCC  = NC  (ID) 

IE(MCTR-NCC) 571, '70,570 

571  NCC  = MCTR 

570  GO  10(550,600,600) ,NSKIP 
600  IE(LSEG-NCC) 561,560,560 

560  L NDSEG  =NCC 
GO  TO  564 

561  L NDSEG  =L  S EG 

564  GO  TO( 555,5673 ,NDTSCO 

567  I E ( NCU  T-L  NDS  rr>)  D66 , 566, 5 65 
566  L NDSEG  =NCUT 

NniSCO=l 

565  CONTINUE 

I F( LIMLOW-NCUT)  311,314,311 

314  DO  315  L P = 1 , L NAX° 

OO  315  NP=1,NNAXP 

315  BS(LP,  NP)  = QS ( LP , N °)  -DELT  AD (L P , NP) 

311  IE(LlMLON-LSEG)5:9 ,568,569 

568  J=2 
NPOINT  =♦ 

L NDSEG  = MCTP 

569  LlML=LIMLOWf J 

DELTA=LAMBOA ( LT ML ' -L AMBDA (LI  ML  ON) 

328  XLL=LAMBDA  (LIMLOW) 

GO  TO( 339,340) , J 
339  NPOINT  =L NDS rG -LT ML OW ♦ 1 

CALCULATE  THE  CONTRIBUTION  TO  INTEG  FRpH  THIS  SEGMENT 
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11 


THIS  C OS'  T POL  SIS T T ON  CONTROLS  HE  COURSE  OF  THE  INTEGRATION  OF  THI 
SEGMENT  II  SCT?  UP  THE  PARTICULAR  COMBINATION  OF  INTEGRATION 
SCHEMES  TO  HE  US^O 

340  IF(9-N POINT) 209, 209, 210 

210  GO  TO(2S  0, 2 02, 20”,  204,205,206,207,208)  ,NPOINT 

202  ASSIGN  260  TO  ISr 
GO  TO  222 

203  ASSIGN  260  TO  1ST 
GO  TO  225 

204  ASSIGN  260  TO  TST 
GO  TO  228 

205  ASSIGN  204  TO  1ST 
GO  TO  222 

206  ASSIGN  204  TO  1ST 
GO  TO  225 

207  ASSIGN  204  TO  1ST 
GO  TO  228 

208  ASSIGN  206  TO  TST 
GO  TO  225 

209  LLSIMP=NP0INT-3 

I F ( ( NP  OI NT/ 2) * 2 -N POINT) 211,212, 21 2 

211  L SI MP=4 
GO  TO  213 

212  LSI MP=7 

ASSIGN  213  TO  1ST 
GO  TO  228 

213  ASSIGN  259  TO  1ST 
GO  TO  228 

259  IF(LSIMP-LLSIM°) ’15,204, 204 
215  ASSIGN  259  TO  1ST 
LSIHP=  LSI MP*2 
GO  TO  225 

END  OF  THE  CONTROL  SECTION  FOR  A SPFCIFIC  SEGMENT 

THIS  SECTION  is  USED  ONLY  AFTER  NDLP24  HAS  PEEN  REACHED 

761  OELT  A=  .08 

XLL=LAMTDA(NGG-2") 

ASSIGN  228  To  IST 
ASSIGN  281  TO  TTRAP 
ASSIGN  264  TO  IFOR 
CALL  PO-KrR<  PO ,XL  AMO, XL G) 

CALL  PO-KER (PLNX  , XLL , XLC) 

DO  76 7 LP=l,LMAXn 
DO  767  NP=1 , NMAX° 

767  9S(LP, NP) =ARRAY(. P,NP,G) 

TRAPEZOIDAL  RULE 

”22  XLM1=XLL 

XLL  = XL  HI  4-DELTA 


161 


oon  ooooo 


! 


IF ( ABS (XLL- XL  G) 0 0 001) 317,213,218 

217  LAST=2 

218  DO  271  LP=1,LM\X3 
PLNKM1  (L P)  = PL  NK  (_  P) 

DO  271  NP=1,NMAX3 

271  DSM 1 (L  D,  NP)  = OS  ( L 3 , NP) 

CALL  PO_KFP(PLNK,YLL,XLG) 

GO  TO  HRAP,  (791 ,783) 

781  G = I NC  ( G) 

DO  782  LP=1,LMAXP 
DO  782  NP=1,NMAXD 

782  BS ( LP > ND) = A PP  AY ( L D » NP , G ) 

GO  TO  250 

783  GO  T0< 220, 221 ) , LAST 

220  RFA0(9)NGEF,  ( (0S(LP,NP)  , LP=1,LMAXP)  ,NP=1, NMAXP) 

GO  TO  250 

221  DO  272  LP=1»LMAXP 

AVOID  TRUNC AT  TON  AT  UPPFR  LIMIT  OF  INTEGRAL 

(NOTICE  THAT  THIS  DOES  NOT  APPLY  WHEN  THE  KLEIN-NISHINA  FORMULA  IS 
BEING  USED) 

I F ( XLG  .GT.  XLNCO  PLNK  (LP)  =PLNKM1  ( 1) 

DO  272  Np=l , UMAX3 

272  BS ( LP, NP) =0  • 

250  H=DELTA»XLG/2 . 

DO  349  L°=l,l.MAXP 
DO  349  NP=1 , NMAX3 
IF(LP-NP) 350, 331, 349 

350  LPNP=L  P*NP 

I F ( (LPNP/2)  * 2 -L  °MP ) 349 , 35  1 , 35 1 

351  I NTEG  ( LP  , NP  ) = INTrG  (LP,NP)  4-H*  (PLNKMl(LP)*OSMl(LP,NP)fXLMl4-PLN<(LP)* 
1BS(LP,NP)/XLL) 

349  CONTINUE 

GO  TO  1ST,  (204,  *>06, 213, 259, 260,  229,342) 

SIMPSON’*S  RULE 

225  XLH2=XLL 

XLM1=XLM2+DELTA 

xll=xlmi+oplta 

IF  ( ABS  (XLL -XLG)-.  0 0 00  1)5  4 2,5  43,543 

542  L AST=2 

543  DO  273  lP=l,LMAX3 
PLNKM2  (L  P)  = PL  M<  ( L p ) 

OO  273  NP=1 , NMAXP 

273  BSM2(LP.  NP)  =BS(L3,  NP) 

REA0(9)NGrr,  ( (PRM1(LP,NP)  ,LP=1,LMAXP)  ,NP=1, NMAXP) 

CALL  POL  <FR ( 3LN<3 1 , XLM1 , XLG) 

CALL  POLKFR(pLNK,YLL,XLG) 

GO  TO( 227,224) , LAST 

'’2  3 READ(9)NGEE,  ( (BS(  LP,NP)  , LP=1,LMAX3)  , NP=1,  NMAXP) 

GO  TO  251 
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00  274  LP=1,LMAX» 

I F ( XLG  .GT.  X L S C ^ > PINK <LP> = PLNKM 1 ( 1 ) 

00  274  NP=1,NMAX° 

274  3S(LP,NP)=0. 

251  H=0ELTA»XLG/7. 

00  346  _P=1,LM AXP 

00  346  NP=1,NMAX° 

IF(LP-NP)  347,  *UP.,  346 

347  LPNP=LP*NP 

IF  ( (LPNP/2)  "?-LD>|P)  346,248,348 

348  I NTEG  ( I.P,  N°  )=  INTrG  iL P , NP ) W ( PL  N<  M2  ( L P)  * BSM2  (LP  , NP)  / XL M2  +4  . * PL  NKMl 
1(LD)*3SH1(L°,  HP)  /XLMi  + PLNK(LP>*9S(LP,NP)/YLL) 

346  CONTINUE 

GO  TO  1ST,  (204, ?06, 213, 259, 260, 228, 342) 

r 

C FOUR  POINT  RULE 

C 

228  XLM3=XLL 

XLN2=XLM3+PrLTA 
XLM 1 = XLN?+0FL  IA 
XLL  = XL  Ml + DELT  A 

1 Ft A8S (XLL-XLG) 0 0001) 54  4,545,545 

544  LAST=2 

545  DO  275  LP  = 1,LUAX° 

PLNKM3  C_P)  =°L  NK  ('.  P) 

DO  275  NP-1,NMAXP 

275  BS93(LP»NP)=9S(L°*NP) 

CALL  pOL  KER ( pLNKu  2 , XL  M2, XLG) 

CALL  P0LKCR(PLNK‘<1,XLM1,XLG) 

CALL  PO.KER  (PLN^,  XLL,  XLG) 

GO  TO  IFOR,  (764  , 269) 

764  G = I NC  ( G) 

S 1 = INC  (3) 

G2=INC  (31 ) 

00  768  Lp=i , LMAX3 
DO  768  N p = 1 , NMA  Xp 
BSM2  (L  P,  NP)  = ARRAY  (LP,NP,G) 

BSMKLP,  NP)  =A?RAY(LP,NP,G1) 

768  3 S ( L 0 , NP)  = APRAv(i_P,NP,G?) 

G = G 2 

I F ( L AS T . FO . 2)  ASSIGN  342  TO  1ST 
GO  T0(252, 277), LAST 

269  READ (9 ) MGEE , ( (9S9 7 (Lp  ,NP)  ,LP=1 ,LMAX3)  ,NP  = 1,MMAXP) 

GO  T 0 ( 257,573)  ,L?  ST 

573  GO  TO  N5K,  ( 266,  257) 

766  NSKIP=N5KIP~1 

ASSIGN  757  TO  NSK 
GO  TO  25  9 

267  REA0(9)NGEE,  ( nS'M  (LP,NP)  ,LP=1,LMAXP>  ,MP=  1,NHAXP) 

GO  T0( 266,574) , LAST 

574  GO  TO(236»27Q),Nr>KID 
1 NSKIP=  NSKIP-1 

GO  TO  257 
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263  GO  T0(  226,2?’')  ,L9S  I 

226  RFAD(9)NGEr,  ( (9S(LP,NP>  ,LP=1,LMXP)  , N P=  1 , NM AXP) 

GO  TO  252 

227  CONTINUE 

00  276  L P=1 , L MAX° 

■ I F ( XLG  .GT.  XLNCO  PL NK  (LP)  = PL  NKM 1 ( 1 ) 

00  276  NP=l,NMA'<n 
276  9S(LP, NP) = 0 . 

252  H=0ELTA*XLG*3./3. 

00  34  3 L P= 1 » L MA  X® 

00  343  NP=l,NMAXn 
IF(LP-NP) 344, 745, 743 

344  LPND  = L P*NP 

t F ( (LPNP/2)  *?-L  011  p)  34  3, 34  5,345 

345  INTEG(LP,Nn)=INTrG(LP,NP)  *•  I- ¥ ( ®L  NKM  3 (L  p)  *BSM3(LP,NP>  f XL  M 3 + 3.  * ®L  NKM2 

1 (LP)  *9  542  (LP,  NP)  /XLM2  + 3.*PLNKM1  (LP)  » BSNl  (LP,NP>  /XLMH-PLNK  (LP)  *9S  (L 

2 P , NP)  / XLL ) 

343  CONTINUE 

GO  TO  1ST,  (204,2  06,213,259,  25  0,  223,342) 

IS  THIS  THE  LAST  SEGMENT  IN  THE  INTEGRAL... 

260  GO  TO( 341,24?) , LAST 
342  MCT  R=  MCT  R +1 
H = H /XL  G 
GO  TO  501 

341  L IML0N=L NOSFG 

RETURN  TO  CALCULATE  THE  NEXT  SEGMENT 
• GO  TO  572 

CALCULATE  MOMENT'  IN  THE  DESIRED  PORTION  OF  THE  L,N  3L  4 N E 

501  LL=0 
NN=  0 
L = 0 
N= 0 

20  NP=NU 
L P=L*1 

651  I F ( NGG-2)  36  E , 3f  6,  7 65 
366  DETER=PC  (LP) ‘LAMBDA  (NCUT) 

DELTAS  C.P,N3)  =SPA  r E ( N , L , S IGO , DETER,  0ELT4B  (L,N)  , DELT  A 9 ( l P ♦ 1 , N ) ,SIGM 
lAT(NCUT)  ) 

365  S IGMAG =S IGM AT  (NGG) 

OTERM=XLG*PO(LP> 

STERM=  OT  FPM  + I NTS ( LP, NP) 

BSC  AT  (LP,NP)  = SpA  3 r ( N , L , S I GO  , STERM  , 9 SC  AT  ( L , N ) , BSCAT  ( L 3 1 , N)  , 5 1 GMAG) 
i/(l.-H*3LN><(LP)  /SI  GMAG) 

BONE(L  P,  NP)  = SPA  Cr  (N,L  ,SIG0,0TcRM,  BONE  (L  , N)  ,OONE  (LPM  , N)  , SIGMA  G ) 
BPNEP ( LP, NP) =BSC( T (LP , NP ) -BONE  (LP, ND) 

26  IF(N-NMAX)  16,  17,17 
16  IF(L) 19, 19, 13 
18  L =L -1 
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N = NU 
50  TO  20 
19  LL  = LLH 

NN=  NN* 1 
L=LL 
N = NN 
GO  TO  20 
C 

17  CONTINUE 

USE  OF  ARRAY  FOR  MOMENTS  STORAGE  REDUCES  PUN  TIME 

I F ( NGG  .LT.  NHL)  GO  TO  765 
DO  762  _P=1,LH\X3 
DO  762  NP=1,NMAXP 
762  ARRAY ( L3 , Np , 5 ) = 3 3 C AT { IP , NP) 

G = INC(G> 

I F (NGG  .LT.  NDLD’4)  GC  TO  765 
ASSIGN  761  TO  ISO 
GO  TO  766 

765  WRITE  (9)  NGG,  ( ( DSC  A T (IP , N P)  , L°=  1 , L M A X ° ) , N P = 1 , NM  A XP) 

766  I F ( NGG  .LT.  NOU7)  GO  TO  754 
NOUT=NOJT+10 

WRITE  (8)  NGG,  ( (90‘JFp  <LP,  NP)  ,LP=1,LMAXP)  ,NP  = i,NMAXP) 

ENERGY =FACT0R/L«M90A  (NGG) 

430  CONTINUE 

754  I F ( NGG  .LT.  NPIGGP)  GC  TO  361 
REWIND  S 

126  DO  438  NG=1, 33 
DO  438  I Y = 1,  5 
DO  438  L P = 1 , 1 0 
DO  780  N P = 1 , 1 9 
700  A ( L D,  N P)  = 0 . 

433  SUMN  ( L P,  I Y,  US ) = P. 

$$$$$  • SSSSS 

BEGIN  PHASr  II 
SSSSS  $$!$$ 

RECONSTRUCTION  0r  T»E  PHOTON  DENSITY  DISTRI3  FROM  113 
MOMENTS  BY  T«E  POLYNOMIAL  METHOD 

624  NGCUT= (NCUT-1 )/10 

¥¥¥¥¥¥¥¥¥¥¥¥¥¥¥¥*.  *¥¥¥¥¥¥¥¥¥¥¥¥¥¥¥¥¥**  ¥¥¥¥•'  ¥*¥¥¥¥¥¥¥***¥¥¥¥%  ******* 

NMIJ  = ANGULAR  VARIABLES  (COS  (THETA)  ) 

NYP=  NUMBER  OF  S^ACE  POINTS  (MEAN  FREE  ° ATHS) 

READ (5, 10 03) ( XMU ( MU) ,MU=1,NMU) 

READ(4,1004) (Y(IY) ,IY=1,NYP) 

I F ( EOF (U)  . N F , 0. ) GO  TO  500 

********  * * ******  5 ******************<-*  ****>-  **************  * ¥ * * * » * * * * 
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00  44  MJ  = 1,NHU 
PMU (1, HU) =1 . 

44  PMU (2, MJ) = XHU  (MU) 

412  I HD  = 2 

NUNIT=8 

UNCOLLIDED  °H0TO')  DENSITY 

630  00  48  ir  = l,N*P 
YIY=Y(IY) 

VIRGIN  (IY.'=07rc>o-  EXP(-YTY)*  (SIG0/YI  Y)  **2  / (CLITE^FOR^TF) 
48  CONTINUE 

0NCE-SCATTrRc0  PHOTON  0FNSITY 

GARB-  0?ER0ASIr-0/(CLITE4F0RPIE  ) 

CALL  FIRST 
NGST=MGCUT 
ASSIGN  531  TO  <ST 
GO  TO  F30 

631  CONTINUE 

DO  710  NG=1,NCCUT 
NGG=1*10*NG 

ENERGY =r ACT OP /LA  MB DA (HGG ) 

710  CONTINUE 

00  414  IY=1,NYP 
DO  414  NG=1,NGMAY 
DO  424  MU  = 1 , N MU 
424  ONFSCT (IY,N3, mu) =0. 

414  SUMN(1 ,IY,NG1 =0. 

OO  415  N P= 1 , N MA  X e 

415  A ( 1 , NP ) " 0 . 

408  LST0P=NST0P/7 

I F ( L END  ,F0.  2)  LSTOP=0 

105  LSTOPP=.STOPM 
NEXT  = ( NS  TOP  + 2 ) / 2 

433  CONTINUE 
436  CONTINUE 

RECONSTRUCTION  0r  THE  LFGFNDRE  EXPANSION  COEFFICIENTS 

00  140  LP=1,IST0RP 
L =L  P-i 
TW0L=L 
XL  = L 

101  NEND  = NF<  T-L 

106  NSTARP=NSUMNO (LP) 

I F( NEND-NST ARP) 1 5 0, 138, 188 

lee  00  621  NP^MST ARP, MEND 

REA D( NUN  IT) NnTGG° , LMAXP, NMAXP 
N=NP-1 


166 


ooo  ooo  ooo  ooo  o o 


CALCULATE  COEFFICIENTS  CF  THE  EXPANSION  POLYNOMIALS 
CALL  W Xp  A N 

43  00  39  N3  = 1 , NO  MAX 

NGG  = 1 + 1 0 * NG 

BRING  MJLTTPL E-SCATTER  MOMENTS  IN  F^OM  OFF  LINE  STORAGE 

410  RF AO(NUNIT)  SICS,  ( f 0 SCAT ( L PP, NPP) , L PP  = 1 , L M A XP)  ,NpP=l,  SHAXp ) 
IF(NGS-NGG)  410,  (-0  9, 4 09 
409  XLAMG=LAMDOA (NGG) 

CALL  BAITA 
DO  39  IY  = 1 , N Y p 
YIY=Y ( IY) 

M SUM=  0 . 

00  38  1=1  , NP 

38  WSIJM  = W SUM  +A  (Lp,  I)  *YIY**  ( I - 1 ) 

BWL>R=^rr  A*WSUM 

39  SUMN(Lp»IY,NG)=SJMN(LP»lY,NG)*RWPR 
REMIND  NUNIT 

621  CONTINUE 

NSUMND  (LP)  -NEN9M 

OUTPUT  SUMN  FOR  USE  IN  CONVERGENCE  ANALYSIS 

DO  119  NGGG= 1 ,NGMAX 
NGG=1+10*NGGG 

119  CONTINUE 
189  CONTINUE 

1 F ( NST  OP-NM  AX  ) 1'-  0 , 12  0,120 

120  CONTINUE 

I F ( L-l >46,46,45 

45  LPJ=2 

USING  RECURSION  RELATION  FOR  LEGENO?E  POLYNOMIALS 

00  47  MJ=1, NMU 
TEMP=PM'J  ( ?,  MU) 

PMU  ( 2 , MJ ) = ALEGRC  ( XL»X  MU  (MU), TEMP,  PM  0(1, MU)  ) 

47  PMU  (1,  MJ)  =TEMP 
GO  TO  117 

46  LPJ=LP 
117  CONTINUE 

DO  742  IY=l,NYp 
YIY=Y( IY) 

ALPHA= 1 ....  CALCULATING  THE  WEIGHT  FUNCTION 
F ACT  = Y IY  * ALPHA 

WT  = SIGOSO*  FACT**  (L  - 2)  * ALPHA*  *2 
WT=WT*  EXP  (-FACT) 

DO  411  NG=1 ,NGMA< 

NGG  = 1 ♦ 10*  NG 

16? 


i 
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WTSUH=Wr»SUwH(LP,  I Y » N G ) 


RECONSTRUCTION  0^  THE  ANGULAR  DEPENDENCE  ...  THE  FIRST  TERM  IN 
THIS  EXPANSION  IS  THE  ALL-ANGLE  DENSITY/  (*  *PIF) 

DO  40  MJ  = 1,N«U 

0NEP(I  Y,  NG,  MU  J-ONEP  (IV,NG,MU)  *■  (TW0L+1  . ) *•'  WT SUM*  PMU  (L  3 J , M U ) /FD  R3 1 E 
40  CONTINUE 

I F ( L .EC.  0)  SJMNdjIY  ,NG)=WTSUM 
411  CONTINUE 

IF(L  .NE.  0)  GO  TO  74? 

TSUM= ( LAMBDA ( 11) -XL  AMO) < SUMN (1 ,IY , 1) /?. 

NG1  = 1 
N G 2 = 2 
NGG  1= 1 1 
NGG  2=21 

740  DLM=LAM^DA(NCG,)-LAMRDA(NGG1) 

rSUM=TSUM+DLM“  (SJMN(1,IY  , NGl)  / L AMRO  A ( NGG  1 ) 4-SUMN  ( 1 , IT , NG  2 ) /L  A M PDA  (N 
1 GG2 ) ) /2. 

I F ( NG?  .EO.  NGMAY)  GO  TO  7A1 

NG1=NG  2 

NG2  = NG  2+ i 

NGG  1=NG3  2 

NGG  2 = N GG  2 «•  1 0 

GO  TO  7k0 

MULTIPLY-SCATTEPED  PHOTON  ENERGY  OENSIT^ 

741  ALENGM  (IY)  = TStM*PACTOR 

742  CONTINUE 
140  CONTINUE 

I F ( NST  OP-NM  AX ) 1 lr  , 1 16 , 1 1 6 

115  NSTOP=N3TOD+NACD 
NST  OpP  =NST0  p+ 1 
GO  TO  409 

116  CONTINUE 

OO  748  I Y = 1 * N YP 

COMBINE  MULTTPLY-SCATTEPED  AND  ONCE-SCA  T TrRED  PHOTON  ENERGY 
DENSITIES  TO  OBT"IN  THE  ENERGY  DENSITY  DUE  TO  ALL  SCATTERED  PHO- 
TONS 

748  ALENG(IV) =ALFNG(IY) fALENGM(IY) 

CALCULATE  UNCOLL*  BED  ENERGY  DENSITY  AND  ENERGY  DENSITY  BUILD-UP 
FACTOR 

00  760  IY=1 * NY° 

VTRGE= VTRGINf IY) ' EZERO 
76  0 ALENGM  (IY)  = (A  LEND  (IY)  *-V2RG£)  / VIRGE 
CONTINUE 

989  WRITE(S,  1103)  (Y  (T  NO)  , INT  = 1 , NY  P) 

WRITE*  2*  110  1)  (Y(IND)  ,INE=i,NYP) 
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1100 

1101 

1102 

1103 

1104 

1105 

r 

1000 

1001 

1002 

1003 

1004 

1005 
1007 
1010 
1011 
1012 


WRIT Ete,  lint.)  ( ALEl'GM(TNO)  ,INO=l,NYP) 

WRI TE ( 2,  1101)  F-’P-P 

WRITE  ( 2,  HOD  (ALCNGM  (INC)  ,IND=1,NYP) 

J0=  1 

DO  145  NG=i,NGMAX 

ngg  = imd*ng 

XL  4 MG  = LA  M DP  A ( NGG) 

ENERGY  =f  ACTOp /XLft  MG 

IE( JO  .FO.  n GO  TO  745 

00  746  IY=1,NYP 

I F ( NG  .LF.  NGCUT)  SUMN(1,IY,NG)  =SUMN(1,IY,NG)  ♦ALA'J3l(IY»NG) 
CONTINUE 
GO  TO  145 
DO  743  IY=1,NYP 

SUHN(1 ,IY,NG) =SM^N(1,IY, NG)*XLAMG** ?/ FACTOR 
CONTINUE 

I F ( JO  .NE.  1)  GO  TO  7^7 
JO=  2 

GO  TO  744 
CONTINUE 

ASSIGN  990  TO  KST 
N G S T = N GM  A X 
DO  124  NG  = 1 , NGST 
NGG=1+10*  NG 

ENERGY -r  4 C^OR /L  A ^ R 0 A ( NGG ) 

DO  124  M'J  = 1,MMU 

ONEP(IY,NG,  MU)=0'JEP(IY,NG,MU)  * FACTOR/ ( ENERGY**  2) 

CONTINUE 

GO  TO  KST,  ( 99  0,  671 ) 

REWIND  5 
CONTINUE 
WRITE(  3,  2013) 

WRITE(6,  2017) 

STOP 

FORMAT  d5 .1) 

FORMAT  (5  F 8 . 4) 

’ FORMAT  (4H  EO? ,^7 . 1 , 12H  KFV,1  MFP=  ,^7,7,7H  METERS) 

; FORMAT ( 3 7 X , 1 6HG U T L 0 UP  FACTORS,//, 10X»llHrNERGY ( KEV)  , 14X , 27  H M 
1REE-PATHS  FROM  SOURCE , / , 21X , 5 ( 5 X , F6 . 3 ) > 

FORMAT  (1  3X,  c6 .1  X ,5  ( 2X,  F9.4  ) ) 

; FORMAT  (1 5H  THIS  RUN  IS  AT,F7.0,7H  METERS) 

1 FORMAT  (2  ( El  4 . 7,  ?y  ) ) 

FORMAT  (7 (12, ?X) ) 

’ FORMAT ( 7 ( El  2 . F» IX)  ) 

! FORMAT  ( 1 2 ( F “ . 2,  1<  ) ) 

. FORMAT  (1 0F7 .3) 

: FORMAT (El  2. 5) 

FORMAT  (12) 

i FORMAT (2 A 6,  i X 0 ,3X,F3.0,3X, El 3. 7) 

FORMAT  (914) 

! FORMAT  (F 6 . 4 , 7 14 ) 
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200  0 FORMAT  ( / , , 8HOM F MFP= , F7 . 3 , 6HMETERS , // , ( 15X, 5 (5X, -7 . 3 ) ) ) 

2001  FORMAT  (///, 3H  L=.I2) 

2002  FORMAT  ( 3X  , F 1 1. , 1 X , 5 ( 1 X , Ell . 4)  > 

2004  FORMAT  (1  AH  CO  S(  TH  F T A ) =MU , 6X,  1 1 H FOR  L AM9DA  = , E12.  5,  9H  ( EN  E RGY= , E 1 2. 5 
1, 5H  KEV)  , //) 

2005  FORMAT  (/,I,,2X,  10F11.4,/,  (5X,10E11.'-0  ) 

2006  FORMAT (2  4X , 3 8 H7 1~  GIN  COMPONENT  Or  T Hr  ° H O T O N OENSIT Y, / , 5 0M  PAVEL  EM 

1GTH(CU)  '1 r AN  - PR  EE  - PATHS  FROM  SOURCE,  //,  ( 15X,  5(  3X  , F7.  7)  , /)  ) 

2000  FORMAT  (2AX,  20HMUL  TIPLY-fCATT  EREO  CO  h = OMEN"  , / , 28  X,  32  H 3 H DT  ON?  =>ER(CM 

1 *»3*KE  VASTER  A 01  AN)  ) 

2009  FORMAT  (//,?‘*Y,?nMT0TAL  PHOTON  DENSITY) 

2010  FORMAT  (//,f^H  WAVELFLGTH  niSTANPE  PROM  SOURCE » / , E OH  ( 

1COMPTON  UNITS)  (IN  MFP"S  OE  VIRGIN  PHOTONS)) 

2011  FORMAT (////, 50H  < ‘*4  DIE*  t *■*  2/  (.Y  +*  L*  EYP  (-Y)  ) ) ^LFGEND^E  EXPANSI0*'  COr 

1 FFICIE  NTS  » / » 1 5X>  ’•  7 H (FOR  L = 0,  MOMENTS  THRU  N = 

2 » 1 2 * 17  H WERE  UTILIZED  ) ) 

2012  FORMAT  (24X,  3AHALL -ANGLE  SCATTERED  PHOTON  DENS1  TY  , /,  2 A X , 2 7 H ( PHOTONS 

1 PER  CM**3  Ppp  K r V ) , / , 5 P H ENERGY (KEV)  MEAM-FRE E-FA T MS  c R 

20M  SOURCE,/, ( 1>X, r (5X,F7 . 3) ,/) ) 

2013  FORMAT  (1  HI) 

2014  FORMAT  (//, COX , 11HFOR  L AMT D A = , El 2 . 5 , 9H  ( EM ERGY= , El  2 . 5 

1 » EH  KEV)  , /) 

2015  FORMAT  (?GX,  At*  HO  A EDUCATION  CE  A PHOTCN  DENSITY  DISTRIBUTION,///,  24  H 

1 SOURCE/  POINT  I ~ riROPIC , /, 1 X, E14.7, 22H  PHOTONS  PEP  SECOND  AT,lX,r 
2t4.7,4H  X rV , / / , GO  H MEDIUM/  INFINITE,  HOMOGENEOUS,  ISOTROPIC,  COUPT 

3 ON  SCATTERED, / , 5X ,43HICA0  STANDARD  a:M0S°HERE  AT  AN  ALTITUDE  OF  ,r 

412. 5, 7H  MFTrRS, /,  9X,  16HF  ARTICLE  D EN  S I T Y = ,E12.5,10H  3 ER  CM«»3,/,9X 
5,  14HC0  NS  T I TUE  NT  S / ,7X,41H  ATOM  DENSITY  ELECTRON  DENSITY 

6) 

1 2C 1 6 FORMAT  <1  5 X , ? A 6,  7 X , E 1 3. 7 , 1 1 X,  E 1 3.  7 > 

2017  FORMAT  (9X,47HTHF  MEAN-FREE  PATH  OF  JNCOLLTDED  PHOTONS  IS,F7.3,?H  M 
1ETERS, //, 

2 32X,5HT0rAL,  / , SOX  , 12Hc'NrR3Y  (KEV), OX, 
3 16H WAVELENGTH  ( C I ) * , 1 IX , 22HCR0SS- SECT  I ON  (CM**-1) ,/> 

201 8 FORMAT (24 X, I 7 , 3X, 3 ( F 1 2. f , 1 2X ) ) 

2019  FORMAT  (S  OH  * 1 ~N  (COMPTON  UNIT)  =COMPTON  WAVELENGTH  OF  THE  ELEC 

1 IRON,  / ,3 pH  EfFRGY  (KEV)  =511. 0 06 /WA V EL f NOT H (CU ) ) 

2020  FORMAT  (//, 35X ,4GHMOMENTS  OF  THE  LEGENDRE  cXPANSION  C D EEF ICI E NTS , // 

1 , 5X , GE  H (THOSE  FOR  N-L  GREATER 

6 T H A N ZERO  OR  ODD  HAVE  NCT  BEEN  CALC JL ATFO) , /,? OX, 10H ( N VS  .)) 

2021  FORMAT  (//,  31H  CO''  MULTIFLY-SCATTERrD  PHOTONS) 

2022  FORMAT {// ,?GH  EO°  ALL  SCATTERED  PHOTONS) 

2023  FORMAT  (24X, 21  HON r - SC AT T EP  COMPONENT , / , 28 X , 3 2HPHOT ON  5 °ER (CM* *3* KEV 
l'STERADI  AN)  ) 

2024  FORMAT  (1  SX,  ? 0 HOMr - SC ATT FR  CUT -OFFS, // , 4 3 H WAVELENGTH  MJ 

1 , //) 

2027  FORMAT  (.35H  ANALYTICAL  ALL-ANGLE  FIRST  SCATTER,//) 

2031  FORMAT  ( P 7 H ALL-ANGLE,  ALL-ENERGY  SC AT T ER- -KEV / C H * * 3 , / , 5 ( 2X , T 

114.7) ) 

2034  FORMAT (24X,43HALL -ANGLE  MULT  I PLY -SC  A T TER  ED  PHOTON  DENSITY,/,  ( 1 5 X , 5 
JL(5Xf  F7.3)  ,/)) 

2035  FORMAT  (///, 44H  ENERGY  BUILD-UP  FACTORS  AT  THE  SPACr  POINTS,//, 14X, 


15  (IX, Ell. 4)  ) 
2040  FORMAT  (16A5) 


ooor?  oooo  oooo 


041  FORMAT  Cl  5>A5  , ///) 
ENO 


SUM ROUTINE  S I G(  F'J rRG Y , ST C-MA) 

COMMON  PMOTC  ? ,40)  , ENFRG  (3,40)  , INCO  (3,  4 3)  , CTA  (7,  17)  , 2\°<  ( * , 1 7)  , PNC  3 
1)  , ON7C  3)  ,NXS,  MEL  , T OLE , T CL  K, 3 1 £ , RZER3 , L‘1A  X , XMUMU, XMUL  * M,  T EMP,  Y I Y , OM 
2UtSIGMA5,EIMTHL,FTGO,ONF,XLNCK,L,MtLP,Mf>,VLAMG,BETA,07.ERO, 

3NTRANS  ,N«U, XL  AMO, G ARB , NV P , NCUT , FACTOR, 3 SC  AT ( 10 , 19)  , 

5SIGMAT  (3  31 ) ,L  AMO"1  A (331)  , ONESCT (5 ,33, 21)  , AL  ANG1  ( 5 , 37)  , AL  E NG  ( 5 ) , XMU  ( 
621)  , Y f 5)  , PMIJ(  ”>,21  ) , VIRGIN  (5) , ALENGMC5)  , A ( 10, 19) 

REAL  I NO  0,  INC  OH 

ARITHMETIC  STATEMENT  FUNCTION  FOR  LOG-LOG  INTERPOLATION 

ALINTP  (X  , XI , Y 1,  X?  , Y 2)  = EXP  (ALOG(Vl)  * AL  OG  (X/Xl)  * AL0G(Y?/Y  1 ) / A'.  OG  ( X 27 
1X1)) 

SIGMA*  0. 

no  58  1=1  , NFL 

CONVF=  ON  ( T)  * 1 . r-'1  A 

I T ( ENE RG  Y -T  01  r)  5'’  ,6  2,63 

62  PHOTO* PHOT ( 1 , 1 ) ONVF 
INCOH*  0. 

GO  TO  69 

63  DO  66  K=? , NXG 

. IT  (AOS (EN ER GY -ENF P G ( I , K) ) -T3L  E) 65, 6F, 64 

64  I F ( EMERGY-f NrRG ( T ,K)>67,65,66 

65  PHOTO* PMOT( I, <) *CONVF 
I NCOH= INCO  ( T , K)  * OONVF 
GO  TO  63 

66  CONTINUE 

67  PHOTO* CO NVF*  ALIN""0  (ENERGY  ,ENERG(I,K-1 ) , PHOT  (I, K-i) , ENERG  (I,K)  , PHOT 

1(1, K)) 

INCOH* CON VF*A LIN" ° ( ENERGY  , ENERG (I , K-i ) , INTO ( I, K-i) , ENERG ( I, <) , INCO 

1(1,  K)) 

68  SIGMA=SIGMA fPHOTO +INCOH 
RETURN 

ENO 


SUBROUTINE  A KERN ( A <ERNL , XLP,XL ) 

COMMON  PHOT  (3,-*0)  , EMERG(3,40)  , INCO  (3, 4 3)  ,rTA(3,  17)  , 0 A PK ( 3,1  7)  , nN(3 
1)  ,0N7( 5)  ,NXS,N~L , T OLE , TOL K, PI E , R7ER3 , L M A X , X MU”U , XMUL A M,  T EMF , Y I Y, DM 
211 , GIGM  AG  , SI  MTHL  , r T GO , ONE  , XLNCK  , L , N , L P , N P , XL  AMG , BETA  , 0 7E  RO , 

3NTRAMS ,NMU, XL  AMO. GARB ,NYP ,NCUT , FACT  OR , BSC  AT (10,  10)  , 

5SIGMAT  ,'33  1)  » L AMR"'  A (371)  , ONTSCT  (5 , 33 , 2 1 ) , AL ANG1  ( 5 , 37 ) , AL F NG( 5 ) , XMU < 
621  ) , Y(6)  ,PMU(  2,  21 ) , VIRGIN  (6)  , ALENGM(5)  , A ( 1 0 , 19) 

ALINTP  (X,  XI  , vi,X'7,  Y2)  =EXP  (ALOG  ( Y1 ) + ALOG  ' X/Xl)  * ALOG  (T2/Y1  ) / A.  OG  (X?/ 
1X1)  ) 
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. 


IF(XL-XLP-2 .0  0 001  ) 95,65,6 
85  I F ( XL  .ST.  XLNCK)  GO  TO  5 

853  R = XLP/XL 
. D=XL-XL° 

DN7K=0  . 

00  854  I = 1 , NFL 

854  DNZK=DNZ<*nN'MT> 

GO  TO  55 

5 IF(A3S (VL-XLP) -1. F-05)6,6,53 

6 A KERNL  =0 . 

RETURN 

53  R=XLP/XL 
0=XL-XLP 

OUMMY=  137.  “FORT  M . *R*  (R+2.+D+D) ) /XL? 

IF<DUMMY-T0LK)S,6  , 55 

55  DN7  K = 0 . 

DO  54  1=1, NEL 

INTERPOLATION  TO  FIND  THE  CORRECTION  FACTOR  (CAPKAY) 
DO  57  K=2,l7 

I F ( A9S (DUMMY- rTA(T ,K) ) -TOLK) 58,58,53 

59  IF  (DUMMY- ETA  (I,  K)  >6  0,56,57 

57  CONTINUE 

60  CAP1=CADK(I  , K -1) 

CAP2  = C A3K (I  ,K) 

I F ( CAP  1-CAP  2) SI," 9,61 

61  CAPKAY=ALINYP(")UMMY,ETA(I,K-1)  ,CAP1,ETA (I,K) ,CAP2> 

• GO  TO  54 

58  CAPKAY=CAPX(  T ,<) 

54  DN7K=DN7K+0N7 (T) ' 0 APKAY 

56  AKERNL=3IE*R7FRO  R + * 2*  (F  * XL/  XL  P-0-9  *Q  * 0)  ♦ PNZ  K 
RETURN 

END 


SUBROUTINE  PQLXE- ( P,XLAMP,XLAMG) 

COMMON  PHOT (7 ,4  0)  , ENERG ( 3,4  0) , INC0(  3,  40)  , rTA (3 , 17) , CA°K ( 3, 17)  , ON ( 7 
1 ) , QNZ ( 3)  ,NXS, N^L,  T OLE , TC L K , P I E , RZERD,  L M A X , X MUMU , XMUL  A M , T EMP, Y I Y » DM 
2D,SIG  MAG, SI  NTHL»C'IGC, ONE, XLNCK, L»N,L°,NP,XXXXX, BETA, D7ERO, 

3NTRANS  ,NMU, XL  AiO, G APB , NY  P , NCUT ,FACT DR , BSC  AT ( 10,  19)  , 

5SIGMAT  (3  31)  ,L  AM'nA  (331)  ,ONE$CT  (5 , 33, 21)  , At  ANG1  (5,37)  , AlENG(5  > , XMU  ( 
621)  , Y(5>  , PM<!(  2,  21  ) , VIRGIN  (5)  , AL  ENGM  ( 5 ) , A ( 10,  19) 

DIMENSION  P ( 2 0 ) 

CALL  AKERN(  °(  1)  , <L  AMP,XI  AMG) 

ARG=1.  +XLAMP-XLAMG 
P ( 2 ) =P  (1)  * ARG 
DO  184  L = 2,  L M AX 
XL  = L 

164  P (L+l) =( <2.*XLM .) * APG*F (L) -XL * P ( L-l ) ) 7 ( XL ♦ 1 . ) 

RETURN 
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SUBROUTINE  ONCE 

COMMON  PHOT  ( 3 , 4 P>  , FNEPG  ( 3 , *,0)  , I NCO  ( 3 , V 3 ) , FT  A (3,  17 ) , OS  °'<  ( 3 , 1 M , ON  ( 3 
1)  , ONZ  ( 3)  ,NXS,  NEL,  T GLE , TOL K , PI E , RZ ETO,  L H A X , XMUMU  , XM'IL  A M , T EMP, Y1 Y,  DM 
2U , S IGM  AS  ,SIMTHL,  ' TGO,ONF  , XLNCK  , L , N,  . P , N e , YL  AMS  , BETA  ,3  ZERO  , 

3NTRANS  ,NM'l,XLAMO,  GARB  , NYP , NC’JT  , FACTOR , 9 SC  A T (10,  19), 

5SIGMAT  (331)  ,L  AMB'iA  (331)  v ONESCT  (5,  33 , 2 1 ) , A l ANGl  ( 5 , 33)  , A l ENG  ( 5 ) , XMU< 
621)  , Y(  5)  ,PMU(  2,  21 ) , VIRGIN  (5)  , ALENGM (5 ) , A ( 10 , 13) 

X M=  XMU  M J 

I F ( X MU  MJ  . EO  . 1.)  XMUMU=0. 99959 
IF(XMIJMJ.LT.XMULAM)  GO  TO  485 
SINTNR=5DRT  (6  33  (1  . -XMUMU**  2)  ) 

ONE  = TEM:”,EXP(-YTv*(SlNTHR+SIGMAS*  ( X MJM'J * S TNTHL -XMUL  A M*  SI  NTHR)  / SI  GO 
1)  /SINTHL) /SINTHR 
X MUMU=  XM 
RETURN 
483  ONE=0. 

RETURN 

END 


SUBROUTINE  FIRST 

COMMON  PHOT  (3,40)  , ENERG(3,40>  , I NC  0 ( 3 , 4 0 ) ,FTA(3,  17),33D<(3,17)  , D.N(3 
1)  , ONZ ( 3)  ,NXS, NEL,TOLt,TGLK,PIE,RZERO,LMAX, XMUMU , XMUL A M , T EMP , YIY, DM 
2U,  SIGMAS,  ST  NTHL,  ST  GO,  ONE,  XLNCK,  L,  N,  L P ,Nf>,  XL  AMG,  SETA,  D ZERO, 

3NTRANS ,NMU, XL  AMD, GARB , NY P , NOUT , FACT  OR, BSC  AT (10 , 19)  , 

5SIGMAT  (331)  , L AM0 "*  A (331)  , ONESCT  (5,33,21)  ,ALANG1(5,33),ALENG(5)  ,XMU( 
621)  ,Y(5>  ,PMU(  2,  21)  , VIRGIN  (5)  , ALENGM (3)  , A (10 ,13) 

REAL  LAMBDA 
OO  488  I Y= 1 , N YD 
YIY=Y(IY) 

ALEN=0. 

NG=1 
NGG  = i 
NGOil 

ASSIGN  731  TO  J ON 
730  XLAMG=LAMPPA ( NGG) 

S IGMAG  =S IGMAT  (NG3  ) 

I F ( NGG  .FB.  1)  XLAMG=XLAMO*. 0001 
I F { NGG  .ED.  NCUT)  XL AMG=LAMBOA ( NC JT) - . 0 01 
CALL  A KERN(  XK  ER»  XL  AMO  , XI  A MG) 

XMULAMni ,-XLAMG+<  L AMO 
SINTHL =5 CRT ( A OS ( 1 . -XMUL7 M** 2)  ) 

TEMP1=GARB*XKER/^INTHL 
TEMP=TEMP1/ YI Y 
ASSIGN  703  TO  ION 
ONETOT  =0 . 

MU=  1 


.. 
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XMU1=1 . 

X MUMU  = 1 . 
f)MU=0. 0001 


C ALCUL  AT  I ON  OF  AN  CLF-OEPfNDENT  ONCE-SCA TTFRED  AND  INTEGRATION 
TO  FIND  All - A MOL  r DENSITY 


701 

705 


707 


700 

703 


706 


CALL  ONCE 
HON=ONE 

X MUMU=  XMUMU  -0  M'J 

I F ( XMU  M J ,lr.  XMJLAM)  GO  TO  I ON , ( 7 0 3 , 7 0 ?) 
CALL  ONCE 

ONETOT  =ONFTOT+OMJ  x (ONE+WON)  72. 

I F ( XMU  MJ  .LE.  0 . '<9901)  DMU~ 0.001 
IF(XMUMJ  ,LF.  0.^9  0 01)  DMU  = 0.01 
I F ( XMU  MJ  .LE.  O.^OOOl)  DM'J=  0.02 
IF ( XMUMJ  .IF.  0.5  0001)  DMU=0.1D 
I F ( NGG  .HE.  MGC)  GO  TO  700 

IF(A9S  (XMU(MIJ)-XUU1)  .GT.  0.0  0 01)  CO  TO  700 
ONESCT  (TY,NG,  MU)  - NON/  (2 . * PIF ) 

I F ( M'J  .FO.  NM l.l)  GO  TO  7C6 

M'J=MU«- 1 

KON=ONE 

XMU1=X  M'J  MU 

GO  TO  701 

XMUMU=  XMUMU  4-0  M'  I 

I F ( DMU  .GT.  0.01)  0MU=  0 .01 

CHEC<=  A9S  (XMU M'.J-Y MULAM) 

I F ( CHE  C<  . L 0.01001)  DMU=0.001 

IF( CHECK  ,LT.  0.00101)  DMU  = C HECK 
XMUMU= XMUMU-DMU 

I F ( OMU  .GE.  CHECK)  ASSIGN  70?  TO  ION 
GO  TO  705 

I F ( NGG  ,NE,  N GO ) GO  TO  498 
XMUMU=  XMU (HU) 

CALL  ONCE 
GO  TO  707 

ALANGl  (IY,NG) =ONrT  OT*  XL A MG** 27  FACT OK 


FIND  ALL-ANCLE  ENERGY  DENSITY  DUE  TO  ONCE-SCATTERED  =>HOTONS 

NG=NG+1 
NGC=NG  C*  1 0 
498  CONTINUE 

GO  TO  JON, ( 73 1,73?,  734,735,736,737) 

731  ASSIGN  73?  TO  JON 
NGG  = NGG«-1 
ONET=OHETOT/XLAHG 
OLM=LAM?OA (NGG) -XL AHG 
GO  TO  730 

7 32  A IE N- A IE N ♦DIM" (ONET+ONETOT/XLAMG) 72. 

I F ( NGG  .LT.  NTRANS)  GO  TO  731 
733  NGG=NGG*1U 


a o a o 


ONETl=  0NET0  7/ XL  A M G 
ASSIGN  73*,  TO  J ON 
GO  TO  730 

734  ASSIGN  735  ?Q  JON 
0NET2=0NET0T/XLA^G 
NGG  = NG  G*- 1 0 

GO  TO  730 

735  ASSIGN  735  TO  JON 

ONET3=ONETOT/XLANG 

NGG=NGGM  0 

DLM=L A M3DA (NGG)-XL AMG 
50  TO  730 

736  ASSIGN  737  ^ JO" 

ALEN=ALENf  (0NET1+3.* (ONET 2 *ONET3 ) fONFTOT/XLAMG) * 3. * DLM/S. 

738  NGG=NGG*2 

0LM=LAM30A ( NGG) -XL AMG 
ONET=ONETOT/XLAHG 
GO  TO  730 

737  AIEN  = AIEN+  ( 0 N ET ♦ ONE T OT/XL AMG ) * OL M/ 2 . 

I F ( NGG  .LT.  NCUT)  GO  TO  738 

ALENG(  IY^ALEN^ACTOR 
488  CONTINUE 
RETURN 

2025  FORMAT (1X,E12.5,3X,F6.3) 

END 


SUBROUTINE  WX  PAN 

COMMON  PHOT  (3,40)  , ENERG( 3,40) , IMCO( 3 , 4 3)  , FT A (3, 17) ,C'D<(3,17),DN(3 
i)  , ON7  ( 3)  ,NXS,NEL«TOLE,TCLK,PIE,RZER0,L1AX,XMUMU,  XM'IL  AM,  T E MF , Y I Y , DM 
2U,SIGMA5,STNTHL,3JG0,0NE,XLNCK,L,N,  „P,NP,XLAMS,  BETA,  0 ZERO, 

3NTRANS ,N MU, XL  AMO, GARR,NYP,NCUT, FACTOR, BSC  AT  (10,19), 

5SIGMAT  (331)  , L AM  =P  A (331)  , ONESCT  ( 5,  33, 7 1 ) , ALANG1  (5,33),AlENG(5),XMU( 
G21)  ,Y(5)  ,PMU(?,?1)  , V I RGI N ( 5)  , A L ENGM  ( 5 ) , A (10,13) 

IF(N)4i,  41,  4 2 

41  A(LP,1)  = 1. 

RETURN 

42  NIP=N 
N=N-1 

TW0NP=2*NI° 

CON=2*  (L  «■  N)  +1 

A (LP,NIP  + 1)  = -A(LP,NIP)/TWONP 
I F ( N)  185,185,  138 
186  DO  37  1=1, N 
I M=  I -1 
X 1=  N- 1 M 
II=NIP-IM 

37  A (LP,II) = ( (CON+XT) *A (LP,II)-A(LP, II-l) ) /TWONP 
185  A (L°,  1 ) = CON'  A (LD,  l)/TWOf.'P 
N = Nf  1 
RETtJRN 

20  05  FORMAT  (/ , 1 3 , 2 X,  10  c 1 1 . 4,  / , (5X,  1 OEi  1.  A ) ) 


oooo  oooo 


SUBROUTINE  BA IT  ft 

COMMON  PMDTfT./jn 


1 ) , ON?  ( 3)  , NXS,  M- L • T OLE , Tl' L K»  PI E,  RZERD » LM  A X , XMUMU  , X M'JL * H » T E MP , y I 
?l),SIGMA3,SINTHL,'TGO,ONF,XLN'CK,L,  N,LPt‘JP,XLAMS,  RETA,r,ER0, 
3MTRANS  * SJ H LI , YLA^o,  G ARB , NYP, NCUT , FACT  OR , BSP  A T ( 10,  19)  , 

5SIGMAT  (7  31)  ,LV1!rf.  ( 331)  » ONESCT  <5 , 33, 21 > , ALANGl < 5 , 331  , AL E NG (5 ) , 
6?1 ) ,Y(?)  , PM  Ht  •»,  21  ) , VIRGIN  (5)  , ALENGM  (5  ) , A ( 10,  19) 

CL  I T E=  ?.  997B2BF4-1  0 
SUM  JN=  0. 

00  36  J3=l , NP 

J=JP-1 

LP? J=L  P+ J*  J 

rERMJ=T0PXAL(L°? I- 1 ) * PSP AT (LP, LP2J) 
rERMJ=TERMJ/TORT^  L ( J*JH*L) 

T ERMJ=  TFRMJ  /TORT'.  L (N-J) 

TERMJ=TERMJ/TORIAL (J) 

J J=  J/2 

IF(  JJ*  2 • EO  . J)  GO  TO  36 

TTOM  !-  _T  CDM  I 


> , DN  ( 3 
V IY, HM 


X M U ( 


II-  ujt  i . tJ  . J ) oU  IU  JO 

TERMJ  = -T  ERM J 

SUM  JN  = SJ  M JM  +T ERMj 

BET  A = Q 7ERG*  T QRI A L ( N)  * SUM  JN  / ( CL  IT  E*X 

RETURN 

ENB 


AMG) 


FUNCTION  SPACE (N,L ,SIGO, DETERM, DEI,  0E2,  SIGH) 

XL  = L 

XN=N 

LP=L+1 

XLP=LP 

LENO=l 

LLND=1 

TERM=0 . 

IF(L-N)550»  F5 1,  rr 1 
LFNO=2 

IF(L) 650,690,691 
LLND=2 

IF(N)  1 3«*  , 13A,  133 

GO  TO ( 65  2 , 69  T ) , LL  NO 

TERM=Xt_*  (XN*-XLP)  ' OE1 

GO  TO  t F32,9B3>  , L r NO 

TFRM=7  ERM+XLPX (XM-XL)*0F2 

HE  TERM  =0ETEPM  +T  E°M*  SIGO/ ( XN* (XL+XL3) ) 

S PACE=  OET  ERM/SI GM 

RETURN 

ENO 


FUNCTION  INC  CO 
K=<  + 1 

I F ( K .GI.  25) 
INC  = K 
RETURN 
END 


function  tori  Aim 

IF(K)  32,  32,  33 
TORIAL'l. 

RETURN 

KTOR=l 

00  34  ir  = i,< 
ktor=ktdr*i t 
torial=<tor 

RETURN 

END 
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